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Abstract 


How  much  information  about  the  shape  of  an  object  can  be  inferred  from  its  image? 
In  particular,  can  the  shape  of  an  object  be  reconstructed  by  measuring  the  light  it 
reflects  from  points  on  its  surface?  These  questions  were  raised  by  Horn  [HO70]  who 
formulated  a  set  of  conditions  such  that  the  image  formation  can  be  described  in  terms 
of  a  first  order  partial  differential  equation,  the  image  irra.dia.nce  equation.  In  general, 
an  image  irradiance  equation  has  infinitely  many  solutions.  Thus  constraints  necessary 
to  find  a  unique  solution  need  to  be  identified. 

First  we  study  the  continuous  image  irradiance  equation.  It  is  demonstrated  when 
and  how  the  knowledge  of  the  position  of  edges  on  a  surface  can  be  used  to  reconstruct 
the  surface.  Furthermore  we  show  how  much  about  the  shape  of  a  surface  can  be 
deduced  from  so  called  singular  points.  At  these  points  the  surface  orientation  is 
uniquely  determined  by  the  measured  brightness. 

Then  we  investigate  images  in  which  certain  types  of  silhouettes,  which  we  call 
b- silhouettes,  can  be  detected.  In  particular  we  answer  the  following  question  in  the 
affirmative:  Is  there  a  set  of  constraints  which  assure  that  if  an  image  irradiance 
equation  has  a  solution,  it  is  unique?  To  this  end  we  postulate  three  constraints  upon 
the  image  irradiance  equation  and  prove  that  they  are  sufficient  to  uniquely  reconstruct 
the  surface  from  its  image.  Furthermore  it  is  shown  that  any  two  of  these  constraints 
are  insufficient  to  assure  a  unique  solution  to  an  image  irradiance  equation.  Examples 
are  given  which  illustrate  the  different  issues. 

Finally,  an  overview  of  known  numerical  methods  for  computing  solutions  to  an 
image  irradiance  equation  are  presented. 


Thesis  Supervisor:  Herthold  K.  P.  Horn 
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Chapter  I 


Motivation 


How  much  information  about  the  shape  of  an  object  can  be  inferred  from  its  image? 
We  are  interested  in  a  special  aspect  of  this  question:  the  reconstruction  problem,  which 
is  to  determine  the  shape  of  an  object  from  measurements  of  the  light  reflected  from  its 
surface.  In  general,  there  are  many  surfaces  which  can  give  rise  to  the  same  image.  So, 
we  will  try  to  identify  and  analyze  constraints  such  that  the  shape  of  a  surface  can  be 
uniquely  reconstructed  from  its  image.  Our  search  for  these  constraints  is  guided  by  the 
following  question:  Are  there  any  “properties”  of  an  object  that  are  also  “properties” 
of  its  image  or  vice  versa?  In  fact,  as  we  shall  see,  such  properties  do  exist. 

Our  work  is  based  on  Horn’s  thesis  [HO70].  He  formulated  a  set  of  conditions  (to  be 
discussed  in  the  next  chapter)  which  lead  to  a  relation  between  the  perceived  brightness 
of  a  small  patch  of  a  surface  and  its  normal  vector.  This  relation,  the  image  11  radiance 
equation,  is  a  first  order  partial  differential  equation  (abbreviated  in  the  following  by 
FOPDE)  and  each  of  its  solutions  determines  the  shape  of  an  object  .  The  problem  of 
finding  solutions  to  the  image  irradiance  equation  is  referred  to  in  the  literature  as  the 
shape  from  shading  problem. 

We  will  take  two  approaches  towards  finding  a  solution  to  the  shape  from  shading 
problem  termed  as  the  local  and  the  global  approach.  By  the  local  approach  we  n.^an 
that  only  a  small  patch  of  an  image  is  used  to  determine  the  shape  of  a  surface.  To 
the  contrary,  in  the  global  approach  we  examine  images  in  which  a  silhouette  can  be 
detected  (here  we  refer  to  the  outline  of  an  image  as  a  silhouette). 

Intuitively,  it  seems  clear  that  by  looking  at  an  image  in  which  a  silhouette  can  be 
identified  we  should  be  able  to  conclude  more  about  the  shape  of  a  surface  whose  image 
we  are  analyzing  than  by  just  looking  at  a  little  patch.  We  will  show  that  from  re  tain 
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images  which  contain  a  silhouette  we  can  uniquely  infer  the  shape  of  the  surface  which 
gives  rise  to  that  image.  Unfortunately  the  global  approach  is  not  always  satisfactory; 
there  are  also  many  images  containing  silhouettes  which  could  be  the  images  of  infinitely 
many  different  surfaces.  There  are  also  infinitely  many  surfaces  which  locally  look  the 
same.  So  we  will  determine  conditions  under  which  the  global  approach  is  better  than 
the  local  approach.  Notwithstanding,  one  can  sometimes  draw  interesting  conclusions 
about  the  shape  of  surfaces  which  give  rise  to  the  same  image  by  just  looking  at  a  small 
patch  of  this  image. 

The  local  approach  is  taken  to  an  extreme  when  we  pose  the  following  question: 
What  can  be  deduced  about  the  shape  of  a  surface  from  so-called  singular  points  of  an 
image  irradiance  equation?  At  these  points  the  surface  normal  to  all  solutions  to  such 
an  equation  is  uniquely  determined  by  the  brightness  there.  We  investigate  the  above 
stated  question  for  a  certain  class  of  image  irradiance  equations,  the  so-called  eikonal 
equations,  which  describe  a  variety  of  physical  phenomena.  For  instance,  experimental 
data  suggest  that  the  flux  of  secondary  electrons  in  a  scanning  electron  microscope  can 
be  described  by  an  eikonal  equation  [LAWH60].  By  using  these  secondary  electrons  to 
modulate  the  appropriate  devices,  an  image  of  a  surface  is  created  by  the  microscope. 
Such  an  image  exhibits  shading  [HO70,  pp. 85-87]  and  therefore  to  determine  the 
shape  of  a  surface  from  its  image  one  effectively  has  to  solve  an  eikonal  equation. 
Studying  eikonal  equations,  we  show  that  the  absolute  value  of  the  Gaussian  curvature 
at  a  singular  point  of  all  surfaces  which  give  rise  to  a  particular  image,  is  the  same. 
Furthermore,  assuming  that  the  surface  is  convex  at  a  singular  point,  we  show  that 
its  shape  can  be  uniquely  determined  in  some  neighborhood  of  such  a  point  from  the 
image  intensities  alone. 

The  other  aspect  of  the  shape  from  shading  problem  which  we  explore  is  its  solution 
when  the  image  contains  a  b-silhouette  (which  is  defined  below).  In  this  case  a  global 
approach  is  taken.  Let  us  first  define  the  bounding  contour  of  a  surface:  a  point  P 
is  on  the  bounding  contour  if  the  line  connecting  the  viewer  and  P  grazes  the  surface 
(i.e.,  if  this  line  lies  in  the  tangent  plane  of  P).  Furthermore  we  assume  that  no  two 
parts  of  a  surface  obscure  each  other,  i.e.,  we  assume  that  the  bounding  contour  is  not 
an  occluding  contour.  The  image  (assuming  orthographic  projection)  of  a  bounding 
contour  will  be  called  the  b-silhouette.  The  surface  normal  at  a  point  on  a  bounding 
contour  is  parallel  to  the  normal  vector  to  the  b-silhouette  and  both  vectors  lie  in  the 
same  plane.  Thus,  some  or  all  of  the  first  order  partial  derivatives  of  the  function 
defining  the  surface  are  infinite  for  points  on  the  bounding  contour  (we  will  say  that 
some  components  of  the  surface  gradient  are  singular  along  a  curve). 

Remark:  In  the  context  of  this  report  the  phrase  a  function  f{x,y)  is  singular  at  a 
point  (xo,2/o)  wiH  always  mean  that: 

Jim  f(x,  y)  =  ±oo.  (1.1.1) 

This  terminology  should  not  be  confused  with  the  notion  of  a  singular  solution  which 
we  will  explain  in  section  A. 4. 
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For  example,  the  bounding  contour  of  a  hemisphere  lying  on  a  plane  B  is  a  circle. 
Consider  a  Lambertian  surface  which  has  the  property  that  each  surface  patch  appears 
equally  bright  from  all  viewing  directions.  If  we  look  at  a  Lambertian  hemisphere  such 
that  the  viewer  and  the  light  source  are  at  the  same  point,  its  b-silhouette  can  be 
determined  from  its  image.  In  this  case  the  image  irradiance  equation  describing  the 
imaging  situation  is  singular  and  all  the  surfaces  which  satisfy  such  an  equation  have 
a  bounding  contour.  (We  extend  the  notion  of  a  singular  function  to  equations.) 

On  the  other  hand,  the  existence  and  the  position  of  a  b-silhouette  cannot  be 
determined  from  a  continuous  image  irradiance  equation.  Thus  the  existence  of  a 
surface  which  satisfies  a  continuous  image  irradiance  equation  and  which  has  a  bounding 
contour  does  not  necessarily  imply  that  all  surfaces  which  satisfy  this  equation  have  a 
bounding  contour. 

In  this  report,  Horn’s  work  on  the  shape  from  shading  problem  is  extended. 
He  studied  primarily  images  of  smooth  surfaces  where  the  imaging  situation  can  be 
described  by  a  continuous  image  irradiance  equation  (which  is  defined  rigorously  in 
section  III.l).  Yet  some  objects  have  edges  and  their  surface  normals  assume  different 
values  depending  upon  which  side  an  edge  is  approached  from  (we  will  say  that  the 
surfaces  normals  are  discontinuous  along  an  edge).  Can  an  edge  be  invisible  in  an 
image?  In  other  words,  is  it  possible  that  a  surface  has  an  edge  without  it  producing  a 
discontinuity  in  the  equation?  As  it  turns  out,  discontinuous  solutions  can  arise  from 
continuous  equations.  In  this  case  initial  conditions  (discussed  in  sections  III. 2,  A.6  and 
A.7)  provide  information  about  the  occurrence  of  a  discontinuity.  Additionally,  we  will 
show  (in  section  III. 2. 3)  under  which  conditions  edges  on  a  surface  can  be  used  as  such 
initial  data. 

We  keep  in  mind  that  the  image  irradiance  equation  describes  a  physical  situa¬ 
tion  and  therefore  the  only  solutions  considered  are  those  corresponding  to  (piecewise) 
smooth  surfaces.  (A  rigorous  definition  of  the  notion  smooth  surface  is  given  in  section 
III.l.)  However  there  are  image  irradiance  equations  for  which  no  solution  correspond¬ 
ing  to  a  smooth  surface  exists.  In  particular,  the  existence  of  such  a  solution  to  a 
singular  image  irradiance  equation  is  not  guaranteed.  We  are  primarily  concerned  with 
the  identification  of  constraints  which  allow  one  to  solve  the  reconstruction  problem 
uniquely. 

In  general,  a  PDE  describes  a  class  of  processes  rather  than  a  particular  one. 
Consider,  for  example,  the  Laplace  equation: 


A/  =  0  (1.1.2) 

where  A  denotes  the  Laplace  operator  and  /  a  scalar  field.  This  PDE  constrains  the 
sources  and  the  curl  of  the  field  /  to  be  zero,  but  an  infinite  number  of  different  fields 
exist  which  fall  into  this  category.  Only  when  some  further  conditions  about  /  are 
specified  can  the  solution  to  the  Laplace  equation  be  restricted  to  a  single  one. 

Similarly,  there  are  an  infinite  number  of  different  surfaces  which  satisfy  a  given 
image  irradiance  equation.  Thus,  as  Horn  [1 1070,  1 1 07 f> ]  has  already  observe1  in 
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general  the  image  irradiance  equation  alone  is  not  sufficient  to  solve  the  reconstruction 
problem  uniquely.  It  remained  an  open  question  whether  there  are  any  imaging  situa¬ 
tions  for  which  every  surface  gives  rise  to  a  different  image.  As  indicated  before,  taking 
the  global  approach  towards  finding  a  solution  to  the  shape  from  shading  problem 
allows  us  to  answer  this  question  affirmatively.  To  this  end  we  will  postulate  three 
constraints  upon  the  image  irradiance  equation  and  prove  that  if  these  constraints  are 
known  to  hold,  the  reconstruction  problem  can  be  solved  uniquely. 


We  now  briefly  describe  the  different  chapters  in  this  thesis. 

Chapter  II  is  a  summary  of  the  issues  involved  in  the  shape  from  shading  problem. 
For  a  more  extensive  discussion  of  this  material  see  [HO70],  [H075]  and  [WOOD78]. 

Chapter  III  gives  an  overview  of  the  different  mathematical  problems  involved  in 
solving  a  FOPDE.  The  mathematical  details  used  in  this  chapter  can  be  found  in 
appendix  I.  In  section  III.l  we  define  two  classes  of  image  irradiance  equations,  the 
continuous  and  the  singular  equations.  Section  III. 2  deals  in  detail  with  the  continuous 
image  irradiance  equation.  In  particular  we  exhibit  in  section  III. 2. 2  that  certain 
continuous  image  irradiance  equations  can  be  transformed  into  singular  equations. 
This  implies  that  all  surfaces  which  satisfy  such  an  image  irradiance  equations  have  a 
bounding  contour.  In  section  III. 2.3  we  examine  how  edges  on  a  surface  can  be  used 
to  reconstruct  the  surface.  To  help  us  understand  the  variety  of  integral  surfaces  of  an 
image  irradiance  equation  we  introduce  in  section  III. 2. 4  some  concepts  from  differential 
geometry  and  in  section  III. 2. 5  gradient  space  as  popularized  by  Mackworth  [MAC73] 
and  Horn  [H077].  Section  III. 3  describes  the  basic  issues  involved  in  solving  a  singular 
image  irradiance  equation.  In  section  III. 3.1  Marr’s  [MA77]  work  on  occluding  contours 
is  reviewed.  Section  III.3.2  discusses  the  method  of  characteristic  curves  for  singular 
image  irradiance  equations. 

In  chapter  IV  we  show  how  a  singular  point  of  an  eikonal  equation  constrains  its 
possible  solutions. 

In  chapter  V  we  postulate  three  constraints  upon  an  image  irradiance  equation  such 
that  the  reconstruction  problem  can  be  solved  uniquely.  In  section  V.l  it  is  proven  that 
if  an  image  irradiance  equation  satisfies  these  constraints  it  has  a  unique  solution.  It 
is  shown  in  section  V.2  that  any  two  of  these  constraints  are  insufficient  to  assure  a 
unique  solution  to  an  image  irradiance  equation. 

Chapter  VI  gives  a  brief  description  and  discussion  of  known  numerical  methods 
for  computing  solutions  to  an  image  irradiance  equation. 

Chapter  VII  summarizes  the  results  of  this  report  and  suggests  some  possible 
applications. 
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Chapter  II 


The  Shape  from  Shading  Problem 


There  are  basically  three  components  to  the  shape  from  shading  problem  which 
have  to  be  taken  into  account.  They  are  the  light  source,  the  object  and  the  camera  as 
depicted  in  figure  1  which  is  taken  from  [WOOD78,  p.32],  and  termed  as  an  imaging 
configuration.  Henceforth  we  will  assume  that  an  image  of  a  surface  is  produced  by 
a  camera.  The  shading  of  such  an  image  can  be  explained  as  follows:  The  exposure 
of  film  in  a  camera  (for  fixed  shutter  speed)  is  proportional  to  image  irradiance,  the 
light  flux  per  unit  area  falling  on  the  image  plane.  Similarly,  grey  levels  measured  in 
an  electronic  imaging  device  are  quantized  measurements  of  image  irradiance.  It  can 
be  shown  that  image  irradiance  in  turn  is  proportional  to  scene  radiance,  the  light  flux 
emitted  by  the  object  per  unit  projected  surface  area  per  unit  solid  angle  [HOS79]. 
The  factor  of  proportionality  depends  on  details  of  the  optical  system,  including  the 
effective  /-number. 

Scene  radiance  depends  on  the 

•  surface  material  and  its  microstructure, 

•  the  incident  light  flux,  and, 

•  the  orientation  of  the  surface. 


Now  we  want  to  relate  the  shape  of  a  surface  to  the  shading  of  its  image,  ~>nsider 
a  viewer-oriented  coordinate  system  with  the  viewer  located  far  above  the  surface  on 
the  2-axis.  If  the  objects  imaged  are  small  compared  to  their  distance  from  the  viewer, 
one  can  approximate  the  imaging  situation  by  an  orthographic  projection: 
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LIGHT 

SOURCE 


Figure  1.  Imaging  configuration 


Zq  Zq 


(II.l.l) 


where  (x,y)  are  the  coordinates  of  the  image  of  a  point  (x,  y,z)  made  with  a  system 
of  effective  focal  length  /,  and  the  viewer  is  at  a  distance  zq  above  the  origin.  We 
assume  that  (x2  -f  y2  +  z2)  <  z$.  For  simplicity  and  without  loss  of  generality  it  is 
also  assumed  that  the  viewing  direction  coincides  with  the  z-axis. 

The  orientation  of  a  patch  of  a  surface  can  be  specified  by  its  gradient  (p,q, — 1), 
where  p  and  q  are  the  first  order  partial  derivatives  of  z  with  respect  to  x  and  y.  For  a 
given  surface  material  and  known  incident  light  flux,  scene  radiance  will  depend  only  on 
surface  gradient.  The  function  which  describes  this  dependence,  R(p,q),  (or  a  contour 
representation  in  gradient  space),  is  called  the  reflectance  map. 

Recall  that  image  irradiance  and  scene  radiance  are  proportional  and  that  we 
assume  orthographic  projection.  If  E(x,y)  is  the  observed  image  irradiance  at  the  point 
(x,y)  in  the  image  then: 

R(p,q)  =  E(x,  y)  (n.1.2) 

where  (p,  q)  are  two  components  of  the  gradient  at  the  corresponding  point  on  the  object 
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being  imaged.  This  equation  is  called  the  image  irradiance  equation.  It  is  clearly  a  first 
order  partial  differential  equation  since  it  involves  only  the  first  order  partial  derivatives 
p  and  q  and  the  coordinates  x  and  y.  In  summary,  to  derive  an  image  irradiance 
equation  we  have  to  know  the  reflectivity  function  and  the  geometry  relating  the  light 
source,  the  object  and  the  camera. 

How  can  an  image  irradiance  equation  be  used  to  analyze  images?  Informally,  sup¬ 
pose  enough  information  about  the  imaging  situation  is  known  so  that  the  reflectance 
map  can  be  derived.  Then  at  every  point  (zo»yo)  the  image  irradiance  denoted  by  Eq 
can  be  measured.  Using  the  image  irradiance  equation  we  can  determine  the  possible 
values  for  p  and  q  such  that: 

R{p,q)  =  E0.  (II.1.3) 

However,  the  values  for  p  and  q  cannot  be  chosen  independently  as  p  and  q  have  to  be 
the  first  order  partial  derivatives  of  a  function  z  =  z[x,  y )  defining  a  smooth  surface.  In 
general  there  are  many  values  for  p  and  q  which  satisfy  equation  (II. 1.3)  and  represent 
the  components  of  a  surface  gradient.  Consequently,  an  image  irradiance  equation  has 
many  solutions  or  equivalently,  for  a  fixed  imaging  configuration  many  surfaces  will 
give  rise  to  the  same  image.  Thus  we  shall  determine  constraints  under  which  the 
solutions  to  an  image  irradiance  equation  are  restricted  to  a  unique  one. 

A  word  of  caution:  Several  issues  such  as  mutual  illumination,  shadows  or  specularity 
are  not  addressed  here. 
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Chapter  III 


Overview 


III .1 .  Classification  of  the  image  Irrad iance  Equation 

In  this  chapter,  we  will  further  investigate  the  shape  from  shading  problem.  Our 
interests  axe  twofold: 

1)  We  wish  to  find  the  solutions  to  an  image  irradiance  equation. 

2)  We  wish  to  determine  constraints  which  guarantee  a  unique  solution  to  an 

image  irradiance  equation. 

The  results  concerning  these  two  issues  are  different  for  continuous  and  singular  image 
irradiance  equations,  both  of  which  are  defined  later  in  this  section.  First  we  make 
some  general  observations  concerning  an  such  equations  and  their  solutions. 

An  image  irradiance  equation  is  a  first  order  partial  differential  equation  in  two 
variables  x  and  y.  We  are  looking  for  a  solution  z  ~  z[x,y)  (also  called  integral  surface) 
which  is  a  function  of  x  and  y  and  defines  a  surface  from  which  light  is  reflected. 
To  be  precise,  the  class  of  functions  which  we  allow  as  solutions  to  the  FOPDE  has 
to  be  specified  and  we  proceed  now  with  some  relevant  definitions.  Using  standard 
nomenclature,  a  function  f(x,y )  is  said  to  be  of  class  Ck  if  it  has  continuous  k- th  order 
partial  derivatives.  The  following  definition  captures  formally  our  geometrical  intuition 
about  a  surface  in  U?3: 
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Definition:  Let  B  C  9?2  be  a  connected  and  closed  region,  and  let  z:B  *-►  SR3  be 
continuous  on  B.  Then  the  single- valued  function  z  =  z(x,y)  defines  a  smooth 
surface  if  z  is  Cl  at  every  interior  point  of  B. 


For  smooth  surfaces 
surface  by: 


the  unit  surface  normal  is  well  defined  for  every  point  on  the 


(1,0,  §*)X(0,  !,§*)■ 
||(l,0,fi)X(0,l,gj)|| 


(III.l.l) 


A  piecewise  smooth  surface  is  defined  similarly.  Unless  otherwise  stated  we  will  always 
assume  that  the  solutions  to  any  given  FOPDE  are  (piecewise)  smooth  surfaces.  We 
will  denote  by  p  and  q  the  first  order  partial  derivatives  of  z  —  z(x,  y)  with  respect  to 
x  and  y. 


We  shall  exploit  the  following  two  facts  about  an  image  irradiance  equation: 

1)  The  equation  does  not  depend  explicitly  on  z. 

2)  The  equation  involves  only  two  functions  R  and  E,  such  that  R  depends 
only  on  p  and  q  and  E  depends  only  on  x  and  y. 

An  immediate  consequence  of  1  is  that  there  are  no  singular  solutions  to  the  equation, 
i.e.,  the  complete  integral  describes  all  possible  solutions.  (For  an  explanation  of  these 
terms  see  section  A.4).  As  for  2,  we  will  study  image  irradiance  equations  which  fall 
into  either  one  of  the  following  two  categories: 

1)  The  functions  R  and  E  are  Cx.  We  say  that  such  an  image  irradiance 
equation  is  continuous.  This  case  is  discussed  in  section  III. 2. 

2)  The  function  R  is  C1 .  The  function  E  is  singular  in  x  and/or  y  but  for  all 
points  (x,y)  at  which  E{x,y)  assumes  finite  value,  it  is  C1.  We  say  that  such 
an  image  irradiance  equation  is  singular.  This  case  is  discussed  in  section  III. 3. 


Since  the  measured  image  irradiance  is  always  finite,  it  seems  at  first  that  singular 
image  irradiance  equations  are  not  of  practical  interest.  However,  for  our  purposes  we 
can  think  about  such  equations  as  useful  mathematical  constructs,  i.e.,  a  singular  image 
irradiance  equation  can  be  obtained  by  transforming  a  continuous  image  irradiance 
equation  appropriately.  By  appropriately  we  mean  that  the  transformation  is  one-to- 
one  and  onto  and  that  the  solutions  to  the  original  equation  and  the  transformed  one 
are  the  same  (section  111. 2. 2).  For  example,  there  is  a  transformation  between  the 
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continuous  image  irradiance  equation  (III. 1.2)  and  the  singular  equation  (III. 1.3): 


P2+q 2 

l+P3  +  <72 

_2  ,  2 
P 


=  *2  +  y2 

*2  +  y2 
i-(*a+ya)‘ 


(m.i.2) 

(III. 1.3) 


We  still  have  to  explain  why  such  a  transformation  is  useful,  i.e.,  why  it  makes 
more  sense  for  us  to  analyze  the  singular  equation  (III. 1.3)  instead  of  the  continuous 
equation  (III. 1.2).  The  reason  is  that  we  can  gain  some  information  about  the  integral 
surfaces  of  a  singular  image  irradiance  equation  without  first  solving  it.  To  this  end  let 
us  define  the  terms  bounding  contour  and  b- silhouette: 


Definition:  A  point  P  lies  on  the  bounding  contour  of  a  surface  defined  by 
z  =  z(x,  y)  if  the  tangent  plane  at  P  is  perpendicular  to  the  x-y  plane.  The 
b-silhouette  is  the  orthographic  projection  of  the  bounding  contour  onto  the 
x-y  plane. 


We  will  show  below  that  each  integral  surface  of  any  singular  image  irradiance  equation 
has  a  bounding  contour,  a  fact  which  is  not  true  for  continuous  image  irradiance 
equations.  In  section  III. 2. 2  we  will  construct  two  integral  surfaces  of  a  continuous 
image  irradiance  equation  such  that  one  of  them  has  a  bounding  contour  whereas  the 
other  does  not.  Why  do  all  integral  surfaces  of  a  singular  image  irradiance  equation 
have  a  bounding  contour?  Recall  (chapter  I)  that  for  points  on  the  bounding  contour 
p  and/or  q  become  infinite.  Now,  an  image  irradiance  equation  is  singular  if  there 
are  points  (x0,yo)  such  that  the  only  values  for  p  and  q  which  satisfy  the  equation  at 
{xo,yo)  are  infinite.  This  explains  the  previously  stated  question.  Furthermore,  as  the 
b-silhouette  is  the  projection  of  the  bounding  contour  onto  the  image  plane,  the  points 
which  constitute  the  b-silhouette  can  be  uniquely  determined  from  a  singular  image 
irradiance  equation  (section  III. 3). 

We  will  not  be  concerned  with  discontinuous  reflectance  maps  since  these  are  rare. 
They  may  occur  when  dealing  with  specularities,  an  issue  not  addressed  in  this  report. 
(It  is  also  assumed  that  the  reflectance  map  is  not  a  constant  since  then  the  image  will 
also  be  constant.) 

III. 2.  The  Continuous  Image  Irradiance  Equation 

In  the  mathematical  literature,  the  most  studied  FOPDE’s  are  continuous  and 
have  continuous  first  order  partial  derivatives.  An  overview  of  known  results  in  this 
area  can  be  found  in  appendix  I  and  we  summarize  them  in  the  next  few  sections 
only  very  briefly.  In  particular  we  explain  why,  in  general,  additional  information  is 
needed  to  restrict  the  solutions  to  a  given  image  irradiance  equation  to  a  single  one. 
Furthermore  we  address  the  problem  of  whether  one  can  deduce  any  properties  of  the 
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integral  surfaces  of  a  given  image  irradiance  equation  by  just  inspecting  the  equation. 
To  this  end  we  investigate  singular  points,  at  which  the  tangent  plane  to  all  integral 
surfaces  of  an  image  irradiance  equation  is  uniquely  defined  by  the  image  intensity 
(section  III. 2. 2).  We  can  then  show  (chapter  IV)  that  for  a  class  of  image  irradiance 
equations  the  absolute  value  of  the  curvature  of  the  integral  surfaces  at  a  singular  point 
is  uniquely  defined  by  the  measured  brightness  there.  In  addition,  we  discuss  a  class 
of  image  irradiance  equations  for  each  member  of  which  we  can  deduce  that  all  of  its 
integral  surfaces  have  a  bounding  contour. 

The  main  tool  for  solving  a  FOPDE  is  the  construction  of  so  called  characteristic 
curves.  The  following  is  a  short  explanation  of  this  notion.  A  more  detailed  exposition 
can  be  found  in  sections  A.2  and  A.3.  Let 

R(p,q)  =  E(x,y)  (III.2.1) 

be  a  continuous  image  irradiance  equation.  What  kind  of  information  about  an  integral 
surface  (defined  by  z  —  z(x,  y ))  of  the  previous  equation  can  be  deduced  from  it?  For  a 
given  point  P  on  z  the  quantities  p  and  q  are  constrained  by  (III. 2.1)  to  lie  on  a  curve. 
Since  p  and  q  determine  the  normal  (p, q, — 1)  at  P,  equation  (III. 2.1)  constrains  the 
feasible  tangent  planes  at  P  to  a  one-parameter  family,  “enveloping  a  conical  surface 
with  P  as  vertex,  called  the  Monge  cone ”  [COHI62b,  p.75].  The  directions  of  the 
generators  of  a  Monge  cone  are  called  characteristic  directions.  Now,  the  characteristic 
curves  are  those  curves  on  an  integral  surface  which  at  every  point  have  as  their  tangent 
direction  a  characteristic  direction.  The  characteristic  curves  can  be  determined  by 
solving  the  characteristic  equations  which  are  five  ordinary  differential  equations  whose 
solutions  depend  on  initial  values: 

t,=F -  dI  =  F’  %-*>+*'  C"-2-2' 

%  =  -(pf.  +  k)  *— (»*  +  #•,). 

Equivalently,  we  can  say  that  for  distinct  initial  values,  different  characteristic  curves 
are  determined.  Since  each  integral  surface  of  a  FOPDE  is  swept  out  by  characteristic 
curves  (which  we  prove  in  section  A. 3),  initial  conditions  (which  are  discussed  further  in 
the  next  sections)  are  necessary  to  restrict  the  solutions  to  an  image  irradiance  equation 
to  a  single  one;  in  the  case  where  an  image  irradiance  equation  is  continuous,  a  unique 
solution  to  the  problem  can  be  found  provided  that  an  appropriate  initial  strip  defined 
by  x  —  x{t),y  =  y{t),z  —  z(t),p  —  p(t)  and  q  —  q{t)  is  known  (section  A. 6  and 
A.7).  By  appropriate  we  mean  that  this  strip  is  not  a  characteristic  strip  and  that  the 
following  condition  holds: 


dx  dy  dy  dx 
ds  dt  ds  dt 


(III. 2. 3) 
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Depending  upon  the  initial  curve,  the  integral  surface  is  either  continuous  or  discon¬ 
tinuous  as  will  be  discussed  in  gi  eater  detail  later  on. 

111.2.1.  The  Linear  Continuous  Image  Irradiance  Equation 

We  first  review  the  case  where  the  image  irradiance  equation  is  a  linear  FOPDE 
eg.,: 

ap~l~  bg  =  E(x,  y )  (III.2.4) 

where  a  and  6  are  constants.  This  equation  is  of  practical  interest  as  its  linear  reflectance 
map  describes  the  reflectivity  properties  of  the  maria  of  the  moon  where  the  constants 
a  and  b  define  the  direction  towards  the  light  source  (i.e.,  the  sun)  [HO70].  Linear 
FOPDE’s  are  special  cases  of  quasi- linear  FOPDE’s  (in  which  a  and  6  are  functions  of 
x,y  and  z).  However,  an  image  irradiance  equation  cannot  be  a  quasi-linear  FOPDE 
unless  it  is  a  linear  FOPDE.  Furthermore  we  can  show  the  following  lemma: 

Lemma:  The  solutions  to  an  image  irradiance  equation  of  the  form: 

f{ap  +  bq)  =  E{x,y)  (III.2.5) 

where  /  is  a  bijection,  f~1(E(x,y))  is  C1  and  where  a  and  b  are  constants, 
can  be  obtained  by  a  simple  coordinate  transformation  from  the  solutions  to 
an  image  irradiance  equation  of  the  form: 

p  =  E{x,y).  (III.2.6) 

Remark:  A  bijection  is  a  one-to-one  and  onto  function.  In  the  above  lemma,  f~~l 
denotes  the  inverse  function  of  f. 

Proof:  Since  the  function  /  is  a  bijection,  equation  (III.2.5)  and  the  transformed 
equation  (III.2.7): 

ap  +  bq  =  f~l(E( x,  y))  (III.2.7) 

have  the  same  solutions.  We  abbreviate  f~1(E(x,y))  by  E{x,y).  To  prove  the  lemma 
we  distinguish  three  cases  depending  upon  the  values  of  a  and  b. 

Case  1)  a  and  6  0 

The  image  irradiance  equation  is  of  the  form: 

ap  +  bq  =  E(  x,  y).  (III.2.8) 

Now  let  z  ~  z(x,  y)  be  a  solution  to: 


P  =  E[x,y) 


(III  .2.9) 


18 


where  x  and  y  are  defined  by: 


x  =  0(1  +  y) 
y  =  b(x~  y). 


(III. 2. 10) 


Then  z(x,y)  =  z(x(x,y),y(x,y))  is  a  solution  to  (III.2.8). 


First  we  express  x  and  y  in  terms  of  x  and  y: 


*  =  -  +  4 

2  a  2b 

_  i _ y_ 

V~2  a  2b' 

Now  we  determine  the  first  order  partial  derivatives  of  z  with  respect  to  x 
terms  of  the  first  order  partial  derivatives  of  z: 

dz  _  dz  dx  dz  dy  _  dz  1  dz  1 

dx  dx  dx  dy  dx  dx  2a  dy  2 a 

dz  dz  dx  ^  dz  dy  dz  1  dz  1 

dy  dx  dy  dy  dy  dx  2b  dy  2b 

dz  dz  dz 

aai  +  bTy  =  T*=E(x’y) 


Thus: 


dz  ,  ,dz 


which  is  the  same  equation  as  (III. 2. 8). 
Case  2)  a  0  and  6  =  0 


The  image  irradiance  equation  is  of  the  form: 


(III.2.11) 


and  y  in 


(III. 2. 12) 


(III.2.13) 


(III.2.14) 


ap  =  E[x,y). 


(II1.2.15) 


As  a  0,  we  can  write  this  equation  equivalently  as: 


E[x,y) 


(111.2.16) 


which  is  of  the  form  (111.2.6). 
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Case  3)  a  =  0  and  6^0 


The  image  irradiance  equation  is  of  the  form: 

bq  =  E{  x,y). 

As  b  ^  0,  we  can  write  this  equation  equivalently  as: 

E{x,y)  ■ 


9  = 


Now  let  z  =  z( x,  y)  be  a  solution  to: 


P  = 


E[x,y) 


where  x  and  y  are  defined  by: 


i  =  y 
y  =  x. 


Then  z[x,y)  ~  z(x{y),y(x))  is  a  solution  to  (III.2.17). 


(III.2.17) 


(III.2.18) 


(III. 2. 19) 


(III. 2.20) 


Now  we  show  that  the  first  order  partial  derivative  of  z  with  respect  to  y  is  the  first 
order  partial  derivative  of  z  with  respect  to  x: 


dz  _  dz  dx _ 

dy  dx  dy 

Thus: 

dz  _  E(x,  y) 
dy  b 

which  is  the  same  as  equation  (III.2.17). 


(III.2.21) 


(III. 2. 22) 


Hence  it  suffices  to  examine  linear  image  irradiance  equations  of  the  form  (III.2.6).  | 

As  previously  mentioned,  a  FOPDE  has,  in  general,  infinitely  many  solutions. 
What  kinds  of  conditions  can  one  impose  such  that  the  solutions  to  a  FOPDE  are 
restricted  to  a  single  one?  A  solution  to  the  Cauchy  problem  which  is  the  problem  of 
constructing  an  integral  surface  passing  through  any  given  curve  C,  provides  us  with 
one  answer  to  this  question  and  is  stated  in  the  following  theorem: 

Theorem:  Let  p  =  E{ x,y)  a  linear  image  irradiance  equation  and  let  C'  be  an 
initial,  continuous  and  non-characteristic  curve.  If  A  (III. 2. 3)  does  not  vanish 
along  C,  then  there  exists  a  unique,  smooth  integral  surface  through  C. 
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Proof:  The  proof  follows  from  the  existence  and  uniqueness  theorem  for  ordinary 
differential  equations  and  can  be  found  in  [COHI62b,  pp.  145-147]  and  in  section  A.6. 

I 


If  A  vanishes  along  C,  then  C  is  either  a  characteristic  curve  and  the  equation  has 
infinitely  many  solutions,  or  the  integral  surface  does  not  have  continuous  derivatives 
along  C.  In  particular  at  the  end  of  section  A. 6  the  following  lemma  is  shown: 

Lemma:  Let  ap  bq  =  E(x,y)  be  a  linear  FOPDE.  Let  C  be  a  non¬ 
characteristic  initial  curve  for  which  A  vanishes  (III. 2. 3).  Then  the  solutions  to 
this  FOPDE  are  cylindrical  surfaces  perpendicular  to  the  x-y  plane. 

Proof:  The  proof  can  be  found  in  section  A.6.  | 


In  the  previous  theorem  C  was  assumed  to  be  a  continuous  curve.  However,  it  can 
be  shown  that  “any  singularities  of  the  initial  data  propagate  in  the  x-y  plane  along 
the  projection  there  of  the  relevant  characteristic  curve”  [GAR6'1,  p.22].  This  is  not 
surprising  since  characteristic  curves  can  be  viewed  as  branch  curves  along  which  two 
integral  surfaces  meet. 

Let  x  =  x(t)  and  y  =  y(t)  denote  the  parametric  representation  of  a  curve  in 
the  x-y  plane.  A  point  P  =  (xo,yo)  is  called  a  double  point,  if  several  values  of  t 
correspond  to  P.  In  the  case  where  an  initial  curve  C  or  its  projection  onto  the  x-y 
plane  has  double  points,  the  integral  surface  containing  C  has  self-intersections  and 
therefore  z  is  not  a  single-valued  function  of  x  and  y. 

The  projection  of  a  characteristic  curve  onto  the  x-y  plane  is  called  a  base  charac¬ 
teristic.  In  the  case  of  a  linear  jjnage  irradiance  equation  the  base  characteristics  are 
unique.  An  integral  surface  of  a  linear  continuous  image  irradiance  equation  can  have 
a  bounding  contour  and  we  can  show  the  following  lemma: 

Lemma:  Let  ap-{-  bq  —  E(x,  y)  be  a  continuous  image  irradiance  equation  and 

let  an  integral  surface  which  has  a  bounding  contour  be  defined  by  z  =  z(x,  y). 

Then  its  b-silhouette  is  a  base  characteristic. 


Proof:  For  points  on  the  bounding  contour  p  and/or  q  are  singular.  Thus  any  initial 
curve  which  has  a  point  in  common  with  the  bounding  contour  is  discontinuous  and 
the  singularity  propagates  along  a  characteristic  curve.  | 

111.2.2.  The  General  Continuous  Image  Irradiance  Equation 

We  discuss  now  general  continuous  image  irradiance  equations.  A  more  detailed 
exposition  of  this  material  can  be  found  in  section  A. 3.  In  deriving  the  system  of 
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characteristic  equations  (III. 2. 2)  it  is  also  assumed  that  an  integral  surface  z  has  con¬ 
tinuous  second  order  partial  derivatives  and  therefore  that  the  equation  py  =  qx  holds. 
In  [PLI54]  it  is  shown  that  an  integral  surface  can  also  be  built  from  characteristic 
strips  in  the  case  when  only  its  first  order  partial  derivatives  are  continuous. 

The  Cauchy  problem  for  a  general  FOPDE  is  formulated  in  the  same  way  as  for 
linear  FOPDE’s  and  is  discussed  in  more  detail  section  A.7.  Let  C  be  an  initial  curve 
specified  in  parametric  form  by  x  =  x(t),y  =  y(t)  and  z  =  z(t).  Then  p(t)  and  q(t) 
along  C  can  be  determined  by  solving  the  two  equations: 

R(p(t),q(t))  =  E(x(t),y(t))  (HI.2.23) 

|=P(0f  +  «K<)|.  (in.2.24) 

As  (II1.2.23)  is,  in  general,  a  nonlinear  equation  in  p  and  q,  several  solutions  may  be 
possible  for  p(t)  and  q(t)  along  C.  To  “avoid  inessential  reference  to  possible  multiple 
valuedness  of  solutions  for  p  and  q  along  C”  [COHI62b,  p.80]  it  is  assumed  that  p  and 
q  are  also  known  as  initial  data.  For  the  following  discussion  we  assume  that  p  and 
q  are  specified  along  C.  In  other  words  the  initial  data  is  an  initial  strip  (which  we 
denote  by  Ci).  Now  if  A  7^  0  (III. 2.3)  along  C  then  a  unique  integral  surface  may  be 
obtained: 

Theorem:  Let  R(p,  q)  =  E(x,  y)  be  a  continuous  image  irradiance  equation. 

Then  for  every  initial  strip  which  is  not  a  characteristic  strip  and  for  which 

A/0,  there  exists  a  unique  integral  surface  through  this  strip. 


Proof:  The  proof  follows  from  the  existence  and  uniqueness  theorem  for  ordinary 
differential  equations  and  can  be  found  in  [COHI62b,  pp.79-82]  and  in  section  A.7. 
It  has  been  shown  [HA28]  that  if  the  initial  data  does  not  have  continuous  second  order 
derivatives,  then  the  solution  does  not  have  continuous  first  order  derivatives.  | 

An  integral  surface  of  an  image  irradiance  equation  is  a  solution  to  a  Cauchy 
problem  only  under  the  assumption  that  R *  +  R2q  7^  0  (section  A.3).  Let  S  be  an 
integral  surface  of  an  image  irradiance  equation  for  which  the  condition  R £  -f-  R*  7^  0 
holds.  Then  there  exists  an  initial  strip  Ci  such  that  the  solution  to  the  characteristic 
equations  with  C\  as  initial  values  is  the  surface  S.  Thus  the  points  of  a  reflectance 
map  at  which  R?  -f-  R7q  =  0  have  to  be  investigated  separately.  Such  points  are  called 
stationary  points  and  are  defined  as: 

Definition:  A  smooth  function  f(x,y )  has  a  stationary  point  at  (i0,j/o)  if 

/x(*o,Vo)  =  0  (III. 2. 25) 

/y(lO>yo)  =  0. 
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We  discuss  now  the  conditions  under  which  stationary  points  constrain  the  integral 
surfaces  of  an  image  irradiance  equation.  First  we  must  define  critical  elements  of  the 
characteristic  equations  of  an  image  irradiance  equation: 

Definition:  Let  R{p,q)  =  E(x,y)  be  an  image  irradiance  equation.  Then  a 
point  denoted  by  {xo,yotPot<}o)  is  a  critical  element  if  (zo,yo)  is  a  stationary 
point  of  E(x,y)  and  ( po,Qo )  is  a  stationary  point  of  R(p,q). 

The  point  (xo,l/diPo><7o)  is  a  critical  point  of  the  image  irradiance  equation  if  it 
is  a  critical  element  and  if  the  values  (zq,  2/o » Po>  <?o)  satisfy  the  image  irradiance 
equation. 

The  point  (io,  Vo,Pq,  qo)  is  a  singular  point  if  it  is  a  a  critical  point  for  which 
the  values  (po,<7o)  are  uniquely  determined  by  the  values  (xo,j/o)- 

A  point  P  is  an  isolated  critical  element  if  in  some  neighborhood  of  it,  it  is  the 
only  stationary  point  of  E(x,  y)  and  R(p,  q).  Isolated  critical  (singular)  points 
can  be  defined  similarly. 


Hence  at  a  singular  point  the  gradient  of  each  integral  surface  can  be  uniquely  deter¬ 
mined  from  the  measured  brightness.  For  example,  the  lines  x  =  0  and  p  =  0  consist 
entirely  of  critical  points  of  the  image  irradiance  equation  p2  =  x2.  Note  that  these 
points  are  not  singular  as  the  value  for  q  cannot  be  uniquely  determined  at  a  point 
(0,  y)  in  the  x-y  plane.  The  point  ( x,y,p,q )  =  (0,0, 0,0),  however,  is  a  singular  point 
of  the  image  irradiance  equation  p2  q2  =  x2  +  y2 .  Thus  the  tangent  plane  to  each 
integral  surface  of  this  equation  at  the  point  (0, 0,  z)  is  parallel  to  the  x-y  plane.  In 
chapter  IV  we  will  show  how  singular  points  constrain  the  integral  surfaces  of  an  image 
irradiance  equation. 

Discussing  the  Cauchy  problem  further,  the  case  A  =  0  for  a  general  FOPDE  is 
analogous  to  the  case  A  =  0  for  a  linear  FOPDE.  If  an  initial  strip  is  a  characteristic 
strip,  then  the  equation  possesses  infinitely  many  solutions.  Again,  the  question  arises 
as  to  whether  it  is  possible  to  specify  an  initial  strip  Cj,  (i.e.,  a  curve  C  and  p  and 
q  along  it),  such  that  A  —  0  and  C\  is  not  a  characteristic  strip.  Such  a  strip  can 
be  specified  but  “then  there  exists  no  integral  surface  which  contains  this  initial  strip 
and  has  continuous  derivatives  up  to  the  second  order  in  its  neighborh'  M”  [COH162b. 
p.83] .  However,  it  might  be  possible  to  construct  a  surface  through  C\.  Then  the  curve 
C  is  a  singular  curve  for  that  surface  [COHI62b,  p.83],  i.e.,  along  this  curve  the  function 
defining  the  surface  does  not  have  continuous  second  order  partial  derivatives. 

Thus  the  solutions  to  a  continuous  image  irradiance  equation  do  not  necessarily 
have  continuous  partial  derivatives  everywhere.  In  particular  p  and  q  can  be  singular  as 
shown  in  the  example  below.  So,  one  integral  surface  of  a  continuous  image  irradiance 
equation  can  have  a  bounding  contour  although  this  does  not  imply  that  eveij  integral 
surface  of  this  equation  has  a  bounding  contour.  For  instance; 

z(x,y)  =  4x1/4  4  f[y) 


(111.2.26) 
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where  /  is  any  continuous  function,  defines  an  integral  surface  of  the  following  image 
ir radiance  equation: 


(p2  -  l)2  _  (1  -  xty 
(P2  +  1)2  (!  +  **)»’ 


(III.2.27) 


Note  that  for  x  —  0,  p,  the  first  order  partial  derivative  of  z  with  respect  to  x  as 
defined  by  (III.2.26),  is  singular.  Thus  this  surface  has  a  bounding  contour  and  its 
image  contains  a  b-silhouette.  However  E(x,y)  and  ^  are  continuous  along  the  curve 

i  =  0. 

On  the  other  hand,  there  exist  integral  surfaces  of  (III.2.27)  which  do  not  have  a 
bounding  contour.  For  example,  the  integral  surface  defined  by: 


*(*.  y )  =  +  ?(y)  (III.2.28) 

where  g  is  any  continuous  function  is  a  solution  to  (III.2.27). 

There  are,  nevertheless,  continuous  image  irradiance  equations  for  which  all  integral 
surfaces  have  a  bounding  contour  as  shown  in  the  next  lemma.  Any  such  equation  can 
always  be  transformed  into  a  singular  image  irradiance  equation  such  that  the  original 
and  the  transformed  equations  have  the  same  solutions: 


Lemma:  Let  R{p,q)  =  E(x,  y)  be  a  continuous  image  irradiance  equation. 
Then  all  of  its  integral  surfaces  have  a  bounding  contour,  if  there  exist  a  C1 
bijection  g,  a  C1  function  /  and  if  the  following  conditions  hold: 

•  The  reflectance  map  can  be  written  as: 

E{p,  q)  =  g{f{p,  q)).  (III.2.29) 

•  There  exist  finite  values  for  x  and  y  denoted  by  io  aud  2/o  such  that: 

Jfiin  g~1{E{x,y))  =  ^oo.  (III.2.30) 

v-*vo 


Proof:  The  solutions  to  the  continuous  image  irradiance  equation  R(p,q)  =  E(x,y ) 
and  the  equation  f{p,q)  =  g~~1{E(x,y))  are  the  same.  Note  that  the  points  (xo.Vo) 
constitute  the  b-silhouette.  Thus  by  virtue  of  the  preceding  discussion  of  bounding 
contours  and  b-silhouettes,  the  lemma  is  self-evident.  | 


In  section  III. 3  and  chapter  V  we  will  discuss  the  solution  to  the  reconstruction 
problem  in  the  case  where  a  b-silhouette  can  be  uniquely  identified  in  an  image. 

The  last  case  of  interest  here  is  where  an  initial  curve  C  degenerates  to  a  point  P 
(sec  also  section  A. 7  and  note  that  the  following  discussion  is  relevant  only  for  FOPDE’s 
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which  are  not  quasi-linear).  Then  “all  characteristic  curves  through  a  fixed  point  P 
of  x,y,z- space  form  an  integral  surface”  [COHI62b,  p.83] .  Such  a  surface  S  has  a 
conical  singularity  at  P  and  the  Monge  cone  is  the  tangent  cone  to  S  at  P.  A  surface 
constructed  in  such  a  manner  is  called  the  integral  conoid  of  the  partial  differential 
equation.  For  a  given  image  irradiance  equation  and  a  given  point  P  there  exists 
exactly  one  integral  conoid. 

In  summary,  the  method  of  characteristic  curves  can  be  used  to  determine  the 
solutions  to  a  FOPDE.  An  initial  curve  which  lies  on  the  integral  surface  must,  in 
general,  be  specified  in  order  to  restrict  the  solutions  to  a  continuous  image  irradiance 
equation  to  a  single  one.  The  critical  points  of  an  image  irradiance  equation  have  to 
be  analyzed  separately. 


1 1 1.2.3.  Edges  and  Vertices 

In  this  section  we  give  examples  of  how  to  formulate  the  Cauchy  problem  such 
that  its  solution  is  an  integral  surface  containing  an  edge.  As  previously  stated,  if  an 
initial  curve  C  is  discontinuous  then  the  surface  normals  to  the  integral  surface  which 
contains  C  are  discontinuous.  We  give  a  very  simple  example  illustrating  this.  Two 
planes  intersect  along  a  straight  line  as  depicted  in  figure  2.  The  image  irradiance 
equation  is  assumed  to  be  linear  and  the  same  for  both  planes: 


p-f  q  =  l- 


(III. 2. 31) 


We  pose  the  following  questions:  Do  the  planes  have  to  be  oriented  in  a  particular 
way  to  assure  that  the  previous  image  irradiance  equation  holds?  Are  there  any  other 
integral  surfaces  which  contain  the  edge  constituted  by  the  intersection  of  two  such 
planes?  The  answers  to  these  questions  depend  upon  additional  information  about  the 
scene  which  is  available.  Referring  to  figure  2,  if  we  specify  the  curve  K  in  three  space, 
it  can  be  used  as  initial  data  to  uniquely  reconstruct  th-  two  planes  from  the  image 
irradiance  equation.  As  the  image  irradiance  equation  is  linear,  the  curve  T  has  to  be 
a  characteristic.  We  observed  earlier  (section  III. 2.1)  that  characteristic  curves  can  be 
viewed  as  branch  curves  at  which  two  different  integral  surfaces  meet.  Through  every 
curve  which  is  not  a  characteristic  there  exists  exactly  one  smooth  surface.  Thus  if 
an  integral  surface  of  (III. 2. 31)  contains  an  invisible  edge,  its  projection  onto  the  r-y 
plane  is  a  line  whose  slope  is  1.  (The  characteristic  curves  ol  a  linear  image  irradiance 
equation  whose  reflectance  map  is  of  the  form  R[p,q)  —  ap  -f-  bq  where  a  and  6  are 
nonzero  constants,  are  straight  lines  whose  slope  is  proportional  to  £.)  Recall  that  the 
linear  reflectance  map  describes,  for  instance,  the  reflectivity  properties  of  the  maria 
of  the  moon  and  that  the  constants  a  and  b  specify  1  lien  the  direction  towards  the  sun 
(section  III. 2.1).  We  conclude  therefore  that  any  ridges  on  the  moon  which  are  not 
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Figure  2.  Two  intersecting  planes 


visible  are  in  the  direction  of  the  sun. 

On  the  other  hand,  if  only  the  edge  T  is  specified,  the  two  surfaces  (which  may 
not  be  planes)  cannot  be  reconstructed  in  a  unique  way  from  a  linear  image  irradiance 
equation.  This  follows  from  the  fact  that  T  has  to  be  a  characteristic  curve.  For 
example,  the  surfaces  defined  by  the  next  two  equations  are  solutions  to  (III.2.31): 


z(x,  y)  =  i  +  (y  ~  x)2  (III.2.32) 

z(x,  y)  =  x  -f-  sin(y  —  x).  (III.2.33) 

At  their  intersection  these  two  surfaces  form  an  edge  which  is  defined  by: 

x  =  y  2  =  1.  (III. 2. 34) 

Furthermore  the  two  planes  (which  are  also  solutions  to  (III. 2. 31))  defined  by: 

x  -  2y  +  z  =  0  (III. 2.35) 

2z  —  3y  +  z  =  0  (III.2.36) 


give  rise  to  the  same  edge. 

If  an  image  irradiance  equation  is  nonlinear,  however,  there  is  only  a  small  number 
of  surfaces  which  can  give  rise  to  a  specific  edge.  In  particular  an  edge  cannot  be  a 
characteristic  curve,  as  two  integral  surfaces  which  meet  along  a  characteristic  curve 
also  have  the  same  tangent  plane  there.  So  let  T  be  an  edge  given  in  parametric  form. 
Using  equations  (III. 2. 23)  and  (III. 2. 24)  we  can  determine  the  possible  values  for  p  and 
q  along  T.  Since  the  image  irradiance  equation  is  assumed  to  be  nonlinear,  several 
solutions  for  p  and  q  arc  possible,  so  many  different  integral  surfaces  can  be  constructed 
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Figure  3.  Vertex 


each  of  which  has  the  edge  embedded  in  it.  Any  two  of  these  surfaces  will  intersect 
at  an  edge.  Hence  in  the  case  of  a  nonlinear  image  irradiance  equation,  only  a  small 
number  of  different  surfaces  can  give  rise  to  a  particular  edge. 

Let  us  now  specify  a  vertex  V  in  space  from  which  three  straight  lines  emanate 
(referred  to  as  a  corner)  as  shown  in  figure  3.  The  equations  of  the  three  lines  are 
given  in  three  space.  Assume,  as  well,  that  the  following  image  irradiance  equation 
holds  in  all  three  regions: 

p2  +  <?2  =  1.  (III.2.37) 

We  ask  how  many  integral  surfaces  exist  which  have  the  same  shading  and  contain  the 
vertex  and  the  three  lines.  (Note  that  these  lines  are  not  visible.)  Clearly,  all  such 
surfaces  will  have  a  singularity  at  the  vertex,  i.e.,  the  function  z  =  z{x,y)  defining  a 
surface  is  not  differentiable  at  V.  A  priori  there  are  two  ways  of  interpreting  the  lines 
leading  from  the  vertex: 

1)  as  curves  lying  on  a  smooth  surface  or 

2)  as  edges. 


We  will  discuss  both  interpretations. 

Assuming  case  1,  we  are  faced  with  the  problem  of  constructing  a  smooth  (except 
at  V)  surface  S  which  has  a  conical  singularity  at  V  and  has  the  three  lines  embedded 
in  it.  Thus  the  vertex  is  a  degenerate  initial  curve  (see  sections  III. 2. 2  and  A. 7)  and  S 
is  the  integral  conoid.  If  such  a  surface  S  exists,  the  image  irradiance  equation  cannot 
be  linear,  the  generators  of  the  integral  conoid  must  be  straight  lines  and  the  vertex 
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cannot  be  a  singular  point  of  the  image  irradiance  equation.  In  section  III. 2.2  it  was 
shown  that  there  cannot  be  another  surface  besides  the  integral  conoid  which  has  a 
conical  singularity  at  V  as  the  integral  conoid  is  uniquely  determined  by  the  image 
irradiance  equation  and  this  vertex.  So  if  the  integral  conoid  with  its  vertex  at  V  can 
be  constructed  such  that  the  three  lines  lie  on  it,  the  surface  S  is  unique.  In  the  case 
where  V"  is  a  singular  point  or  the  generators  of  the  integral  conoid  are  not  straight 
lines,  no  surface  5  can  be  found  and  the  three  lines  necessarily  specify  three  edges. 

We  show  now  that  to  be  able  to  construct  the  surface  S  as  described  in  the  preceding 
paragraph,  the  three  lines  necessarily  are  characteristic  curves  of  the  image  irradiance 
equation.  So  suppose  that  the  integral  conoid  can  be  constructed  at  V  such  that  it 
contains  the  three  lines.  Then  the  integral  conoid  and  the  Monge  cone  are  identical, 
which  follows  from  the  fact  that  the  Monge  cone  is  the  tangent  cone  to  S  at  V.  Thus 
if  the  integral  conoid  has  straight  lines,  they  must  coincide  with  the  generators  of  the 
Monge  cone.  The  lines  leading  from  the  vertex  are  therefore  characteristic  curves  of 
the  image  irradiance  equation.  In  other  words,  if  the  lines  are  characteristic  curves  and 
if  V  is  not  a  singular  point,  then  an  integral  conoid  which  contains  the  corner  exists. 

In  our  particular  example  the  integral  conoid  of  (III. 2. 37)  is  a  right  circular  cone 
whose  generators  are  inclined  to  the  x-y  plane  at  45  degrees.  (A  cone  is  called  a  right 
circular  cone  if  the  angle  between  its  axis  and  a  plane  containing  the  circular  cross 
section  is  90  degrees.)  Thus  if  the  three  lines  are  characteristic  curves,  the  right  circular 
cone  is  the  desired  surface  S. 

The  other  interpretation  of  the  corner  is  that  the  three  lines  originating  at  the 
vertex  specify  three  edges.  Again  only  in  the  case  where  the  image  irradiance  equation 
is  nonlinear  can  we  possibly  find  three  surfaces  which  form  the  three  edges  since,  if 
the  image  irradiance  equation  is  linear,  different  integral  surfaces  intersect  only  along 
characteristic  curves.  Yet  the  base  characteristics  are  all  parallel  for  a  linear  equation, 
hence  the  projection  of  the  three  lines  onto  the  x-y  plane  cannot  be  base  characteristics 
and  so  the  lines  cannot  be  characteristics. 

We  want  then  to  find  three  surfaces  which  intersect  along  the  three  lines  using  these 
lines  as  initial  data.  For  a  particular  corner,  three  such  surfaces  do  not  necessarily  exist 
as  can  be  easily  seen  using  equation  (III. 2. 37).  The  integral  surfaces  of  (III. 2. 37)  that 
are  constructed  using  a  straight  line  (which  is  not  a  characteristic)  as  initial  data,  are 
planes  inclined  at  45  degrees  to  the  x-y  plane.  Although  through  three  lines  which 
constitute  a  corner  we  can  always  find  three  planes,  they  are  not  necessarily  inclined 
at  45  degrees  to  the  x-y  plane  and  are  therefore  not  necessarily  integral  surfaces  of 
(III. 2. 37).  However,  if  three  integral  surfaces  do  exist,  they  are  unique.  From  the  two 
lines  which  are  embedded  in  every  surface,  p  and  q  can  be  uniquely  determined.  If 
these  values  for  p  and  q  as  functions  of  x  and  y  satisfy  the  image  irradiance  equation 
then  each  of  the  surfaces  can  be  uniquely  determined. 

Summarizing  the  observations  made  above:  if  three  lines  which  form  a  corner  are 
characteristic  curves  of  an  image  irradiance  equation,  then  the  surface  containing  the 
corner  is  the  integral  conoid  of  the  equation.  Otherwise  such  a  corner  can  be  used  to 


reconstruct  the  three  surfaces  which  form  it. 

111.2.4.  Gaussian  Image 

In  the  previous  sections  we  reviewed  the  well  known  results  concerning  continuous 
FOPDE’s  and  showed  that  in  general,  a  continuous  FOPDE  has  infinitely  many  solu¬ 
tions.  Only  after  imposing  some  additional  constraints  (in  particular  by  specifying  an 
initial  strip  which  is  not  a  characteristic  strip)  can  the  solutions  be  restricted  to  a  single 
one. 

It  is  clear,  then,  that  a  continuous  image  irradiance  equation  does  not  contain 
sufficient  information  to  uniquely  reconstruct  the  shape  of  a  surface  from  its  image. 
What  partial  knowledge  about  a  surface  does  an  image  irradiance  equation  provide?  To 
understand  this  problem  better,  certain  concepts  taken  from  differential  geometry  are 
now  introduced  as  they  provide  us  with  a  formalism  to  discuss  the  shape  of  a  surface. 
Some  technical  prerequisites  are  first  reviewed  briefly.  A  more  detailed  exposition  can 
be  found  in  any  standard  book  on  differential  geometry,  e.g.,  [CAR76].  For  the  rest  of 
this  section  we  assume  that  a  function  z  —  z(x,  y)  which  defines  a  surface  is  C 2. 

The  Gaussian  sphere  is  a  sphere  of  unit  radius.  The  surface  normal  at  every  point 
on  the  Gaussian  sphere  can  have  two  possible  orientations,  referred  to  as  the  positive  (or 
outward)  and  negative  (or  inward)  directions.  We  adopt  the  convention  that  a  surface 
normal  on  the  Gaussian  sphere  points  outward.  The  Gaussian  mapping  maps  a  point  P 
on  a  surface  into  a  corresponding  point  P°  on  the  Gaussian  sphere  such  that  P  and  PG 
have  the  same  surface  normal.  This  mapping  is  well  defined  since  no  two  points  on  the 
Gaussian  sphere  have  the  same  normal.  The  Gaussian  image  (sometimes  referred  to  as 
the  spherical  image)  of  a  connected  patch  of  a  smooth  surface  maps  into  a  connected 
region  on  the  Gaussian  sphere.  The  area  of  the  Gaussian  image  of  a  surface  is  called 
its  integral  curvature. 

Points  on  a  surface  are  classified  according  to  the  behavior  of  the  tangent  planes 
with  respect  to  the  surface: 

Definition:  A  point  P  on  a  surface  is  elliptic  if  the  tangent  plane  at  P  does  not 
intersect  the  surface  at  any  other  point.  A  point  P  is  hyperbolic  if  the  tangent 
plane  at  P  intersects  the  surface  in  a  curve  which  has  two  branches  intersecting 
at  P.  A  point  P  is  parabolic  if  the  tangent  plane  at  P  intersects  the  surface 
along  a  single  curve  (which  can  degenerate  to  a  point). 

A  surface  which  consists  only  of  elliptic  points  is  locally  convex.  We  will  say  a 
surface  is  hyperbolic  if  it  consists  only  of  hyperbolic  points.  The  points  which  separate 
regions  where  a  surface  consists  of  elliptic  and  hyperbolic  points,  respectively,  are 
parabolic  points. 

The  grouping  of  points  on  a  surface  into  elliptic,  hyperbolic  and  parabolic  points 
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Figure  4.  Gaussian  image  of  a  neighborhood  of  an  elliptic  point 


can  easily  be  expressed  in  terms  of  Gaussian  curvature.  Let  F  be  a  simply  connected 
and  bounded  surface  patch  which  is  enclosed  by  the  dosed  curve  k  and  let  G  denote 
the  Gaussian  image  of  F,  which  is  enclosed  by  the  closed  curve  kG.  “We  divide  the 
area  G  enclosed  by  kG  on  the  sphere  by  the  area  F  enclosed  by  k  on  the  surface  and 
then  shrink  the  curve  k  down  to  a  point  P  on  the  surface.  Then  F  and  G  approach 
zero  and  their  quotient  approaches  a  definite  limit  K: 

lim  ~=K.  (III.2.38) 

F 

The  number  K  defined  in  this  way  is  called  the  Gaussian  curvature”  [HICV52,  pp.193- 
194].  A  point  P  on  a  surface  is  elliptic  if  the  Gaussian  curvature  is  positive  there, 
hyperbolic  if  the  Gaussian  curvature  is  negative,  and  parabolic  if  the  Gaussian  curvature 
is  zero.  If  z  =  z(x,y)  specifies  3  surface  and  is  C2,  then  the  Gaussian  curvature  at  a 
point  ( x ,  y,  z)  is  defined  by: 


K{x,y) 


ZxzZyy  %Xy 

(1  +  zl  +  xfp’ 


(III. 2.39) 


The  common  notions  of  a  surface  being  convex  or  concave  do  not  refer  to  the  type 
(as  defined  above)  of  points  a  surface  consists  of.  They  merely  distinguish  the  two  sides 
of  a  locally  convex  surface  (or  correspond  to  the  two  possible  directions  of  a  normal 
vector)  with  respect  to  a  viewer. 

The  mapping  between  a  surface  and  the  Gaussian  sphere  is  one-to-one  if  the  surface 
is  either  locally  convex  or  hyperbolic,  “If  we  move  around  an  elliptic  point  along  a 
small  closed  curve  that  lies  on  the  surface,  its  spherical  image  -  assuming  that  the 
surface  has  no  double  points  -  will  also  be  a  closed  curve  without  double  points,  and 
this  curve  is  traversed  in  the  same  sense  as  the  original  curve  (figure  4).  A  small  curve 
without  double  points  about  a  hyperbolic  point  is  also  mapped  into  a  curve  without 


30 


Figure  5.  Gaussian  image  of  a  neighborhood  of  a  hyperbolic  point 


double  points,  but  in  this  case  the  sense  is  reversed  (figure  5)”  [HICV52,  pp. 195-196]. 
The  spherical  image  of  a  surface  which  consists  of  elliptic,  hyperbolic  and  parabolic 
points  consists  of  several  sheets,  i.e.,  several  points  on  a  surface  get  mapped  into  the 
same  point  on  the  Gaussian  sphere. 

We  will  also  need  the  definitions  of  a  closed  and  a  compact  surface: 

Definition  [CAR76,  p.l  12):  Let  A  be  a  subset  of  5ft3.  We  say  that  p  6  5R3  is  a 
limit  point  of  A  if  every  neighborhood  of  p  in  5R3  contains  a  point  of  A  distinct 
from  p.  A  is  said  to  be  closed  if  it  contains  all  of  its  limit  points.  A  is  bounded 
if  it  is  contained  in  some  ball  of  5R3.  If  A  is  closed  and  bounded  it  is  called  a 
compact  set. 

For  example,  the  surface  of  a  sphere  is  compact,  whereas  a  paraboloid  of  revolution 
defined  by  z(x,  y)  =  x2  -f-  y2  is  a  closed  but  not  compact  surface. 

The  notion  similar  captures  formally  what  we  mean  by  saying  that  two  surfaces 
have  the  same  shape: 

Definition:  Two  surfaces  in  5R3  are  similar  if  they  can  be  mapped  into  each 
other  by  a  composition  of  translations,  rotations,  reflections  and  dilations. 

III. 2.5.  Gradient  Space 

In  this  section  we  will  briefly  discuss  gradient  space  as  popularized  by  Mackworth 
[MAC73]  and  Horn  [H077]  and  which  is  now  a  standard  tool  in  vision  research.  Our 
objective  is  to  show  how  the  concepts  from  differential  geometry  help  us  to  understand 
the  variety  of  integral  surfaces  of  an  image  irradiance  equation.  We  will  also  investigate 
how  different  constraints  may  restrict  the  number  of  possible  solutions  to  a  given  image 
irradiance  equation. 
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Figure  6.  Map:  reflectance  map  -*  Gaussian  sphere 

Let 

R(p,  q)  —  E{x,  y)  (III.2.40) 

be  an  image  irradiance  equation.  For  simplicity  of  exposition  we  assume  that  E{ x,y) 
is  defined  at  every  point  in  the  x-y  plane.  Then  we  can  write  the  previous  equation  as 
a  system  of  two  equations: 

R{p ,  q)  —  c  (III.2.41) 

E[x ,  y)  —  c  (III.2.42) 

where  c  is  a  constant.  In  the  p-q  plane  (also  called  gradient  space),  the  graph  of 
R(p,q)  =  c,  for  all  possible  values  of  c,  is  called  the  reflectance  map.  A  reflectance 
map  can  be  mapped  onto  the  Gaussian  sphere  by  placing  the  south  pole  of  the  Gaussian 
sphere  onto  the  origin  of  the  p-q  plane  by  means  of  a  simple  projection  from  the  center 
of  the  sphere,  as  illustrated  in  figure  6,  Note  that  the  mapping  between  gradient  space 
and  the  Gaussian  sphere  is  conformal.  The  graph  of  E(x,  y)  —  c,  for  all  possible 
values  of  c,  can  be  drawn  in  the  x-y  plane,  referred  to  as  the  image  plane.  For  example 
consider  the  following  eikonal  equation: 

p2  -j-  q2  =  x2  4  y2.  (III. 2. 43) 

This  equation  can  be  rewritten  as: 

p2  -f  q2  =  c 

S  +  y2 


=  c. 


(III. 2. 44) 
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Figure  7.  Graphs  for  the  equation  p2  -f-  q2  —  x2  -f-  y 2 


The  graphs  for  these  equations  are  depicted  in  figure  7. 

Let  z  —  z(x,  y)  define  an  integral  surface  of  an  image  irradiance  equation.  Then 
the  gradient  at  every  point  of  z  is  the  vector  ( p ,  q,  — l)  where  p  =  ||  and  q  =  |^.  As 
discussed  in  chapter  II,  we  assume  that  each  point  P  —  (x0,  Vo,  z{x0,  y0 ))  on  the  surface 
is  mapped  via  orthographic  projection  into  the  point  ( x0 ,  y0 )  in  the  image  plane.  We  also 
define  a  mapping  from  surface  orientation  to  gradient  space:  The  gradient  (p0)q0, — 1) 
at  P  is  mapped  into  the  point  (po,qo)  in  gradient  space. 

An  image  irradiance  equation  gives  us  some  information  about  the  correspondence 
between  points  in  the  image  plane  and  points  in  gradient  space.  In  particular,  if  (io,  2/o) 
lies  on  the  curve  E(x,y)  =  c o  for  some  constant  co,  then  (po,<7o)  lies  on  the  curve 
q)  =  cq-  Note  that  determining  an  integral  surface  (up  to  a  constant  factor)  of 
an  image  irradiance  equation  is  equivalent  to  specifying  the  correspondence  between 
points  in  the  image  plane  and  points  in  gradient  space. 

We  want  to  show  now  that  any  two  integral  surfaces  of  an  image  irradiance  equation 
need  not  be  similar.  For  two  such  surfaces  to  be  similar,  it  is  necessary  that  they  consist 
of  the  same  type  (i.e.,  elliptic,  hyperbolic  or  parabolic)  of  points.  An  image  irradiance 
equation,  however,  does  not  restrict  the  type  of  points  on  the  integral  surfaces.  So 
by  knowing  one  integral  surface  of  an  image  irradiance  equation,  other  such  surfaces 
cannot,  in  general,  be  obtained  through  similarity  transformations.  To  prove  these 
assertions  we  examine  some  integral  surfaces  of  (111.2.43).  In  particular  the  surface 
defined  by  the  following  equation  consists  entirely  of  elliptic  points: 


z{x,y)  = 


x2  -f  y2 


2 


(III. 2. 45) 


Figure  8.  Graph  for  characteristic  curves 


whereas  the  surface  defined  by: 


is  hyperbolic. 


z  ~  xy 


(III. 2.46)  , 


As  already  noted  several  times,  the  basic  tool  for  solving  a  FOPDE  is  the  method 
of  characteristic  curves.  Using  gradient  space,  the  characteristic  curves  of  an  image 
irradiance  equation  and  the  necessity  of  specifying  an  initial  curve  to  restrict  the  possible 
solutions  to  a  single  one  can  be  visualized  easily  [H075],  In  the  following  we  use 
only  four  of  the  five  characteristic  equations  (III.2.2)  of  an  image  irradiance  equation 

r{p,  9)  =  Eix>  y): 

dJL  =  RJr>.a\  (III. 2.47) 


—  =  flP(P.9) 

~  =  Ex(  x,y) 
as 

j-  =  Ey{x,y). 
as 


(111. 2. 48) 

(111. 2. 49) 

(111.2.50) 
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In  the  image  plane,  dx  and  dy  can  be  viewed  as  the  two  components  of  a  vector  and 
Ex  and  Ey  as  the  direction  of  a  normal  vector  to  the  curve  E(x,y )  =  c  (where  c  is 
a  constant).  In  the  same  way  dp  and  dq  can  be  interpreted  in  gradient  space  as  the 
components  of  a  vector  and  Rp  and  Rq  as  the  direction  of  a  normal  vector  to  the  iso¬ 
brightness  curves  R(p,  q)  =  c.  From  equations  (III.2.47)  and  (III. 2. 48)  one  can  deduce 
that  the  two  vectors  (dx,dy)  and  [Rp,Rq)  are  parallel.  Similarly,  equations  (III.2.49) 
and  (III.2.50)  indicate  that  the  two  vectors  ( dp,dq )  and  ( Ex,Ey )  are  parallel. 

Now  let  the  point  (20,2/0)  in  the  image  plane  correspond  to  the  point  (po,9o)  in 
gradient  space  (see  also  figure  8  which  is  taken  from  [WOOD78,  p.l90j).  Then  a 
step  in  the  image  plane  from  the  point  (20,2/0)  in  the  direction  of  a  characteristic 
curve,  corresponds  to  a  step  in  gradient  space  from  (po,<?o)  in  the  direction  ( Ex,Ey ). 
Conversely,  a  step  from  the  point  (po,<7o)  in  the  direction  ( dp,dq )  corresponds  to  a 
movement  in  the  image  plane  starting  at  (20,2/0).  in  the  direction  (Rp,Rq).  Specifying 
an  initial  curve  now  gives  a  correspondence  between  a  set  of  points,  denoted  by  X ,  in 
the  x-y  plane  and  a  set  of  points,  denoted  by  P,  in  gradient  space.  If  this  initial  curve 
is  not  a  characteristic,  the  characteristic  curves  can  be  expanded  from  each  point  in  X 
which  just  corresponds  to  expanding  curves  in  gradient  space  from  every  point  in  P.  In 
this  manner  a  point  (p,  q)  is  assigned  to  each  point  in  the  image  plane.  Once  p  and  q  are 
known  at  every  point  on  the  surface,  the  function  z  —  z{ x,y)  defining  the  surface  can 
be  found  by  integration  (assuming  that  p  and  q  satisfy  the  equation  py  =  qx).  Thus  to 
determine  a  unique  (up  to  translation  in  the  2-direction)  solution  to  an  image  irradiance 
equation,  enough  information  has  to  be  known  so  that  the  characteristic  curves  can  be 
expanded  simultaneously  in  the  image  plane  and  in  gradient  space. 


III. 3.  The  Singular  Image  Irradiance  Equation 

In  section  III.  1  we  classified  image  irradiance  equations  into  two  categories.  In  the 
previous  sections,  we  discussed  equations  which  fall  into  the  first  category,  i.e.,  the 
case  where  E(z,y)  is  C1.  This  section  deals  with  the  second  class,  i.e.,  the  case  where 
E(x,y)  is  singular. 

Recall  (chapter  I)  that  a  function  E(x,y)  is  singular  at  a  point  (20,2/0)  if-’ 

Jjrn  E{x,y)  =  ±00.  (III.3.1) 

v-*vo 


Lemma:  Let  R(p,q)  =  E{ x,y)  be  a  singular  image  irradiance  equation.  Then 
the  b-silhouette  consists  of  those  points  (20,2/0)  in  the  x-y  plane  for  which 
E{x,y)  is  singular. 


Proof:  Each  singular  image  irradiance  equation  defines  a  b-silhouette.  The  lemma 
follows.  | 


Let  us  recall  here  briefly  the  geometry  of  the  image  forming  system  as  discussed 
in  chapter  II.  It  is  assumed  that  the  viewing  direction  coincides  with  the  z- axis. 
Furthermore,  we  approximate  the  imaging  situation  by  orthographic  projection  hence 
the  lines  connecting  the  viewer  and  points  on  the  surface  are  parallel  to  each  other  and 
perpendicular  to  the  x-y  plane.  For  simplicity  we  will  restrict  our  attention  to  images 
where  a  single  b-silhouette  occurs.  Images  with  multiple  b-silhouettes  can  be  treated 
in  a  similar  fashion.  We  will  call  a  b-silhouette  nondegenerate  if  it  is  a  connected  and 
smooth  curve  in  the  x-y  plane. 

Let  R[p,q)  =  E[x,y)  be  a  singular  image  irradiance  equation,  let  the  b-silhouette 
be  defined  by  w(x,y)  =  0  and  let  z  =  z(x,y)  define  an  integral  surface.  Then 
its  bounding  contour  consists  of  the  points  (x0,  yo,  2(2:0,  yo))  on  the  integral  surface 
such  that  (z0,t/0)  belongs  to  the  b-silhouette.  We  assume  that  integral  surfaces  are 
(piecewise)  smooth  (section  III.  1).  So,  no  parts  of  an  integral  surface  obscure  each 
other  and  therefore  the  bounding  contour  is  a  set  of  (piecewise)  continuous  curves.  As 
previously  stated,  the  lines  connecting  the  viewer  and  points  on  the  bounding  contour 
graze  the  surface.  In  other  words,  at  every  point  on  the  bounding  contour  the  tangent 
plane  is  perpendicular  to  the  x-y  plane.  In  terms  of  p  and  q  this  means  that  p  and/or  q 
assume  infinite  value  at  points  on  the  bounding  contour.  The  spherical  image  of  points 
on  the  bounding  contour  (section  III. 2.4)  corresponds  to  points  on  the  equator  of  the 
Gaussian  sphere. 

If  the  equation  of  the  b-silhouette  is  known,  the  surface  normal  to  points  on  the 
bounding  contour  can  be  determined.  (This  follows  directly  from  the  definition  of 
bounding  contour.  The  surface  normal  at  a  point  (x0,  y0,  z(x0,  yo))  on  the  bounding 
contour  is  (±u'x(zo,  yo),  iu,v(Io,  Vo),  0)  where  the  +  and  —  sign  distinguish  the  two 
sides  of  a  surface  with  respect  to  the  viewer.)  This  surface  is  tangent  to  a  cylinder 
whose  intersection  with  the  x-y  plane  is  the  b-silhouette  and  whose  generators  are 
perpendicular  to  the  x-y  plane. 

So  again  we  pose  the  question:  What  are  the  constraints  necessary  to  obtain  a 
unique  solution  to  a  singular  image  irradiance  equation?  We  shall  proceed  in  two  ways. 
In  section  III. 3. 2  we  show  that  the  method  of  characteristic  curves  can  be  used  to  find 
the  integral  surfaces  of  a  singular  image  irradiance  equation.  Unfortunately,  specifying 
a  bounding  contour  does  not  restrict  the  possible  integral  surfaces  to  a  single  one,  as 
shown  by  means  of  an  example  in  section  III. 3. 2.  In  chapter  V  we  investigate  constraints 
which  enable  us  to  solve  the  reconstruction  problem  uniquely. 

The  existence  of  a  smooth  integral  surface  is  not  guaranteed  for  every  singular 
image  irradiance  equation  as  we  now'  demonstrate.  It  is  important  to  notice  that  we 
are  looking  for  a  global  solution  to  an  image  irradiance  equation,  i.e. ,  an  integral  surface 
which  is  defined  at  every  point  ( x,y )  for  which  the  equation  is  defined.  In  particular, 
an  integral  surface  should  be  bounded  for  points  (x,y)  which  lie  on  the  b-silbouette.  An 
example  of  an  image  irradiance  equation  for  which  no  global  bounded  solution  exists 
is  given  by: 
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p  =  -.  (111.3.2) 

x 

The  general  solution  to  this  equation  is: 

z(x,y)  =  lnx-f  /(y)-fc  (ill.3.3) 

where  /  is  any  continuous  function  and  c  is  a  constant.  For  points  on  the  y-axis  (i.e., 
x  =  0)  z(x,y)  as  defined  by  the  previous  equation  is  not  bounded. 


For  any  given  singular  image  irradiance  equation,  the  points  on  the  b-silbouette 
can  be  found  by  inspection  of  E(x,y).  The  following  example  should  clarify  this: 


2  i  2 

P  +9  — 


x2  4-  y2 
i  —  (x2  4-  y2)' 


(III.3.4) 


Along  the  circle  x2  4"  2/2  —  1,  £(x,y)  assumes  infinite  value.  An  integral  surface  of  the 
previous  equation  is  a  sphere  defined  by: 


z(x,  y)  =  y^l  _  (x2  +  y2)  (111.3.5) 

with  p  and  q  given  by: 


p  =  — - - -  (III.3.6) 

V1  —  (*2  +  y2) 

- . 

\/l  —  (x2  +  y2) 

Along  the  circle  x2  +  y2  =  1,  p  and  q  are  infinite  but  the  surface  normal  there  is  well 
defined  as  (x,  y,  0).  Note  also  that  if  z  =  z(x,  y )  is  a  solution  to  an  image  irradiance 
equation,  then  so  is  z  =  z(x,y)  -f-  c  for  any  constant  c.  As  we  are  solely  interested  in 
determining  the  shape  of  surfaces,  we  will  assume  in  the  following  chapters  that  c  is 
given.  (In  other  words,  all  integral  surfaces  are  determined  only  up  to  translation  in 
the  z-direction.  Thus  when  we  say,  for  instance,  that  there  is  a  unique  integral  surface 
of  an  image  irradiance  equation  we  mean  that  it  is  unique  up  to  translation  in  the 
z-direction.)  In  section  V.2  we  show  that  the  convex  and  the  concave  hemispheres  are 
the  only  two  integral  surfaces  of  (III. 3. 4)  which  have  continuous  second  order  partial 
derivatives. 

111.3.1.  Marr’s  Occluding  Contours 

One  of  the  first  to  investigate  b-silhouettes  was  David  Marr  [MA77].  His  nomencla¬ 
ture  differs  from  ours  and  to  avoid  confusion  we  compare  the  two  terminologies.  In  his 


37 


Figure  9.  Generalized  cone 


paper  bounding  contours  are  referred  to  as  contour  generators,  and  b-silhouettes,  as 
contours. 

Marr  poses  the  question  of  what  can  be  inferred  about  the  shape  of  a  surface 
from  its  b-silhouette  alone.  In  order  to  more  easily  approach  this  problem,  a  priori 
assumptions  about  the  surface  must  be  made.  The  restrictions  imposed  are: 

Rl.)  The  surface  is  defined  by  a  C 2  function. 

R2.)  Each  point  on  the  bounding  contour  projects  to  a  different  point  on  the 

b-silhouette. 

R3.)  Nearby  points  on  the  b-silhouette  arise  from  nearby  points  on  the  bound¬ 
ing  contour. 

R4.)  The  bounding  contour  is  planar. 

As  discussed  in  the  paper,  the  first  three  restrictions  are  not  severe.  Marr  points  out 
that  for  bounded  surfaces,  R3  follows  from  R2  but  he  nevertheless  chooses  to  state 
the  two  restrictions  separately  as  “they  have  sufficiently  different  meaning”  [MA77]. 
The  third  restriction  states  that  there  are  no  gaps  in  the  viewing  direction.  Implied 
by  this  is  that  a  bounding  contour  cannot  be  created  by  two  surfaces,  one  of  which 
partially  occludes  the  other,  as  this  would  violate  R2.  A  consequence  of  the  first  three 
restrictions  is  that  a  bounding  contour  is  a  continuous  curve. 

The  main  technique  used  in  [MA77]  to  interpret  b-silhouettes  is  to  examine  their 
inflexion  points,  i.e.,  those  points  which  separate  convex  regions  on  a  b-silhouette  from 
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concave  ones.  If  a  bounding  contour  is  assumed  to  be  planar,  inflexion  points  on  the 
b-silhouette  imply  the  existence  of  inflexion  points  on  the  bounding  contour.  This 
observation  leads  to  restriction  R4  which  as  Marr  observes,  is  a  very  strong  one.  He  is 
therefore  only  able  to  analyze  a  small  number  of  scenes. 

The  main  theorem  of  [MA77]  (theorem  1)  states  that  generalized  cones  are  the  only 
type  of  surfaces  which  satisfy  the  four  restrictions  above  plus  an  additional  assumption 
about  the  viewing  direction.  “A  generalized  cone  (as  shown  in  figure  9)  is  the  surface 
generated  by  moving  a  smooth  cross  section  p  along  a  straight  axis  A.  The  cross  section 
may  vary  smoothly  in  size  (as  prescribed  by  the  axial  scaling  function  h(z)),  but  its 
shape  remains  constant”  [MA77].  Not  only  is  R4  a  very  strong  assumption,  but  for 
theorem  1  to  hold,  the  viewing  direction  is  confined  to  a  plane  which  lies  parallel  to  a 
cross  section  of  the  cone  whose  shape  we  wish  to  determine.  However,  a  priori,  there 
is  no  way  to  determine  if  the  viewing  direction  satisfies  the  constraint  of  the  theorem! 

The  proof  of  theorem  2  as  stated  in  [MA77]  is  wrong.  This  theorem  claims  that 
a  surface  obeys  the  four  previously  stated  restrictions  for  all  viewing  directions  if  and 
only  if  it  is  a  quadratic  surface.  In  the  proof  given  it  is  assumed  that  the  surface  can 
be  defined  by  a  polynomial;  the  theorem  is  therefore  only  shown  to  hold  for  this  special 
case. 

1 1 1.3.2.  The  Method  of  Characteristic  Curves 

Let  us  recall  here  that  the  method  of  characteristic  curves  can  be  used  to  solve 
a  continuous  image  irradiance  equation.  Is  this  method  still  valid  in  the  case  where 
an  equation  is  singular?  In  fact,  the  answer  is  positive.  The  method  of  characteristic 
curves  is  based  on  a  local  theory  as  stated  in  [COHI62b,  p.62] :  “It  should  be  emphasized 
again  that  all  statements  and  derivations  are  in  the  small,  i.e.,  they  concern  merely 
neighborhoods  of  points,  etc.,  without  necessarily  specifying  the  extension  of  these 
neighborhoods.”  It  is  precisely  the  local  nature  of  this  theory  which  allows  us  to  apply 
it  to  singular  image  irradiance  equations. 

Let  R(p,  q)  —  E(x,  y)  be  a  fixed,  singular  image  irradiance  equation.  The  charac¬ 
teristic  curves  can  be  constructed  in  a  fashion  analogous  to  the  continuous  case.  A 
difficulty  arises  only  when  the  Cauchy  problem  is  stated  for  this  equation.  In  the  con¬ 
tinuous  case,  for  any  strip  denoted  by  C\  (as  defined  in  section  III. 2. 2)  which  is  not 
a  characteristic  strip  and  for  which  A  7^  0  (III .3.3),  an  integral  surface  can  be  found 
which  has  C\  embedded  in  it.  Furthermore  this  integral  surface  is  uniquely  specified  by 
the  FOPDE  and  the  strip  C\.  The  theorem  used  to  obtain  this  uniqueness  result  is  the 
existence  and  uniqueness  theorem  for  ordinary  differential  equations.  Unfortunately, 
in  some  neighborhood  of  a  b-silhouette  this  theorem  does  not  hold  anymore,  so  there 
is  no  guarantee  that  the  equation  will  have  a  solution  for  any  initial  conditions. 

To  understand  why  the  existence  and  uniqueness  theorem  for  ordinary  differential 
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Figure  10.  Solutions  to  m.3.8 


equations  does  not  hold  at  a  singularity,  we  examine  the  equation: 

dy  _  x-j-y 
dx  x 

which  is  singular  at  (x,y)  —  (0,0).  The  set  of  solutions  are: 

y  =  x  ln(ci)  (III.3.8) 

where  c  denotes  any  constant.  The  graphs  of  some  of  the  curves  defined  by  the  previous 
equation  are  depicted  in  figure  10  [SMI68a,  p .  1 48] . 

The  initial  condition  that  the  solution  must  pass  through  the  point  P,  which  has 
the  coordinates  (0,0),  does  not  suffice  to  pin  down  the  solution  to  (III.3.7)  uniquely. 
Furthermore  there  is  no  solution  which  passes  through  the  point  P°  with  the  coordinates 
(0,3).  This  implies  that  for  the  initial  condition  P°  no  solution  exists.  In  summary, 
one  cannot  claim  that  for  any  initial  condition,  (III.3.7)  can  be  solved  in  a  unique  way. 
Moreover,  the  type  of  singularity  in  an  ordinary  differential  equation  constrains  what 
initial  data  is  consistent. 

If  an  image  irradiance  equation  is  singular,  i.e. ,  if  E(x,y )  is  singular,  then  some 
or  all  of  its  characteristic  equations  are  singular.  Thus  when  specifying  a  bounding 
contour,  there  is  no  guarantee  that  a  solution  exists  which  has  the  given  bounding 
contour  embedded  in  it.  Furthermore  we  can  show  the  following  lemma: 

Lemma:  A  singular  image  irradiance  equation  can  have  several  integral  surfaces 
each  of  which  has  the  same  bounding  contour. 


(III.3.7) 
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Proof:  We  will  prove  the  lemma  by  showing  it  for  a  particular  image  irradiance  equation. 
The  two  surfaces  defined  by  (II1.3.10)  are  solutions  to  (111.3.9)  and  each  has  the  same 
bounding  contour: 


Pi  +  q1=  4  +  1  (ni.3.9) 

*t 

z(x,  y)  =  —  4i*  4-  y  (III.3.10) 

z{x,y)=  4z*-fy. 


I 


In  chapter  V  we  shall  prove  that  the  only  integral  surfaces  of  (III. 3. 4)  are  the 
convex  and  concave  hemispheres.  So  for  any  bounding  contour  C  other  than  the  one 
specified  by  i2  -f-y2  =  1  ,z  =  0,  no  integral  surface  exists  which  has  C  embedded  in  it. 
Note  in  addition  that  the  concave  and  the  convex  hemispheres  have  the  same  bounding 
contour.  By  assuming  that  the  surface  is  convex  in  the  direction  of  the  viewer,  the 
solutions  to  (III.3.4)  are  restricted  to  a  single  one  although  this  is  not  true  if  we  forego 
this  assumption,  as  shown  in  the  proof  of  the  previous  lemma. 
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Chapter  IV 


Singular  Points 


IV.  1.  Basic  Concepts 

The  question  which  we  analyze  in  this  chapter  is:  How  much  information  about 
the  shape  of  a  surface  can  one  obtain  from  a  singular  point  of  an  eikonal  equation?  At 
a  singular  point  of  an  image  irradiance  equation  the  gradient  of  all  its  integral  surfaces 
is  defined  by  the  image  intensities  there  (section  III. 2. 2).  As  discussed  in  chapter  I, 
equations  of  the  form: 

p2  -f  q2  =  E(x,y)  (IV.1.1) 

are  called  eikonal  equations  and  describe  for  instance  the  flux  of  the  secondary  electrons 
in  a  scanning  electron  microscope  since  this  varies  approximately  as  /(p2  -f-  q2)  where 
/  is  a  continuous  function  [LAWH60]. 

To  obtain  our  results  we  will  impose  some  technical  conditions  upon  E(x,  y)  (which 
we  will  discuss  later  in  this  section)  and  shall  refer  to  eikonal  equations  which  satisfy 
these  conditions  as  constrained  eikonal  equations.  The  two  results  which  we  prove  are: 

•  There  exist  exactly  two  locally  convex  integral  surfaces  of  a  constrained 

eikonal  equation  in  some  neighborhood  of  a  singular  point. 

•  At  a  singular  point,  the  Gaussian  curvature  of  each  integral  surface  of  a 

constrained  eikonal  equation  has  the  same  absolute  value. 

The  first  statement  can  be  expressed  in  other  words  as:  if  z  =  z[x,  y)  defines  one 
locally  convex  solution,  then  z  —  —  z(z,y)  defines  the  other.  Hence  this  result  can  be 
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viewed  as  a  uniqueness  result  modulo  the  concave/convex  ambiguity.  The  second  result 
exhibits  a  situation  where  a  “property”  of  surfaces  can  be  determined  from  their  given 
(common)  image.  Our  theorems  apply  also  to  more  general  equations:  in  appendix  II 
we  exhibit  a  class  of  image  irradiance  equations  whose  solutions  can  be  obtained  by 
solving  an  appropriate  eikonal  equation. 

To  characterize  a  constrained  eikonal  equation  we  need  the  following  definition: 

Definition:  A  function  g[x,y)  vanishes  precisely  to  second  order  at  the  point 
(0,0)  if  its  (limited)  Taylor  series  expansion  at  (0,0)  is  of  the  form: 

g{x,  y)  =  ctx2  +  f3xy  -f  nfy2  -f  o((|xj  -f  jy|)2)  (IV.1.2) 

where  a,/?  and  q  are  constants,  at  least  one  of  which  is  nonzero. 


Then  a  constrained  eikonal  equation  is  defined  as: 

Definition:  An  eikonal  equation  p2  +  q2  =  E(x,y)  is  constrained  if  E[x,y)  is 
a  C3  function  satisfying  the  following  conditions  in  some  neighborhood  of  the 
point  (x0,  yo)- 

1)  (io,yo)  is  a  stationary  point  of  E(x,y) 

2)  E(x0,y  o)=0  (IV.  1.3) 

3)  E[x,  y)  >  0  for  (x,y)  ^  (x0,yQ) 

4)  E(x,y)  vanishes  precisely  to  second  order  at  (xo.yo)- 


Let  us  discuss  these  conditions  a  bit  further.  The  reflectance  map  of  an  eikonal 
equation  is  R{p,q)  —  p2  -f-  q2  and  therefore  its  stationary  point  is  given  by  p  —  0  and 
q  —  0.  We  have  imposed  the  condition  that  (xo,yo)  he  a  stationary  point  of  E(x,  y). 
Thus  the  point  P  —  {x,y,p,q)  —  (xo>2A)>0, 0)  »s  a  critical  point  of  a  constrained  eikonal 
equation.  Whence  it  follows  from  condition  2  that  P  is  a  singular  point  of  such  an 
equation.  From  the  third  condition  we  can  deduce  now  that  P  is  an  isolated  singular 
point  (section  III.2.2).  By  using  a  suitable  linear  transformation  we  may  assume, 
without  loss  of  generality,  that  the  point  (0,0)  is  the  stationary  point  of  E{x,y).  Since 
E(x,  y)  is  assumed  to  be  positive  near  the  origin: 

ax2  +  0xy  +  7y2  >  0  for  [x,  y)  ^  (0, 0)  (IV.  1.4) 

defines  a  positive  bilinear  form.  Thus  the  subsequent  inequality  [BRSET5,  p.  1 82)  holds: 

B 2 

c*7 - —  >0. 


(IV. 1.5) 
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Moreover,  the  constants  a  and  7  in  the  (limited)  Taylor  series  expansion  of  E(x,y)  must 
be  positive. 

As  mentioned  earlier,  if  z  =  z( x,  y)  defines  a  locally  convex  solution,  then  so  does 
z  =  —z(x,y).  Hence,  to  simplify  subsequent  discussions,  we  will  say  that  z(x,y )  is 
a  locally  convex  solution  to  a  constrained  eikonal  equation  if  it  satisfies  the  following 
positivity  conditions  in  some  neighborhood  of  the  origin: 

1)  *(0,0)  =  0 

2)  z(x,y)€C3  (IV.1.6) 

3)  z(x,y)  >  o. 

Throughout  this  chapter  it  is  assumed  that  a  solution  z  ==  z(x,  y)  satisfies  the  positivity 
conditions  and  that  the  origin  is  the  isolated  singular  point  of  a  constrained  eikonal 
equation. 

Horn  [H075]  used  singular  points  to  compute  initial  conditions  sufficient  to  solve 
an  image  irradiance  equation.  In  particular  he  assumed  that  an  integral  surface  is 
convex  at  a  singular  point.  In  general,  however,  such  a  surface  does  not  exist.  Our 
results  show  when  Horn’s  method  can  be  used  and  we  prove  that  in  those  cases,  no 
initial  conditions  are  needed  to  compute  the  convex  surface. 

IV. 2.  Preliminaries 

The  results  we  prove  in  this  chapter  require  a  number  of  technical  prerequisites 
which  are  introduced  in  this  section.  One  of  the  key  concepts  is  that  of  a  Taylor  series, 
which  is  discussed  in  order  to  deal  with  the  problem  of  approximating  a  function  by 
polynomials.  To  be  able  to  write  the  Taylor  series  expansion  of  a  Ck  function  z(x,  y)  in 
a  concise  form  we  introduce  the.notion  of  a  homogeneous  polynomial,  which  is  a  sum 
of  terms  of  the  same  degree: 

Definition:  A  polynomial  P(x,  y)  is  a  homogeneous  polynomial  of  degree  k  if  it 

is  of  the  form: 

k 

P(x,  y)=^MV-y  (IV.2.1) 

>=0 

where  for  each  j,  hj  denotes  a  constant. 

Definition:  Let  z(x,y)  be  a  Ck  function.  Then  its  (limited)  Taylor  series 

expansion  at  the  origin  is: 

k 

z(x>  V)  =  zo  +  Zj  +  o((|x|  +  |t/|)fc) 


(IV.2.2) 
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where  zq  denotes  a  constant  and  for  each  j,  Zj  is  a  homogeneous  polynomial  of 
degree  j.  For  any  j  such  that  0  <  j  <  k,  there  exists  exactly  one  polynomial 
Zj  which  matches  z  at  the  origin  up  to  the  j-th  derivative. 


Let  us  denote  by  /  the  first  k  -(-  1  terms  in  the  Taylor  series  expansion  of  z(x,y)  at 
the  origin.  Then  /  approximates  z  up  to  fc-th  order,  i.e.,  the  error  between  /  and  z 
vanishes  faster  than  any  polynomial  of  degree  k. 


Now  we  investigate  the  characteristic  equations  of  a  constrained  eikonal  equation 
and  their  solutions  in  some  neighborhood  of  an  isolated  singular  point.  The  four  relevant 
characteristic  equations  are: 


dp 

dt 


Ex(x,y) 


dq 

dt 


Ey(x,y). 


(IV.2.3) 


Since  E{x,y)  vanishes  precisely  to  second  order  at  (0,0),  we  can  rewrite  the  charac¬ 
teristic  equations  in  some  neighborhood  of  the  origin  as: 


dx 

dt 

dy 

dt 

dp 

dt 

dq 

dt 


2  P 
2  q 

2ax  -j-  py  -f  o(|z|  +  |y|) 
fix  +  2^y  +  o(|x|  -j-  )yj). 


(IV.2.4) 


One  can  view  x,  y,p  and  q  as  coordinates  in  9i4.  So,  let  £  denote  the  four- tuple  (x,  y,  p,  q) 
and  let  F(x,  y,  p,  q)  =  p2  -|-  q2  —  E(x,  y).  We  will  say  that  £  £  F  if  x,  y,  p  and  q  satisfy 
the  equation  F  =  0.  Using  this  notation  we  introduce: 


Definition:  Let  p2  -{-  q2  =  E{x,y)  be  a  constrained  eikonal  equation  and  let  £ 
denote  the  four-tuple  (x,  y,p,  q).  The  characteristic  equations  can  be  written 
as: 


d{=M  +  G(i) 


(IV.2.5) 


(IV.2.6) 


(IV.2.7) 


(IV.2.8) 


To  describe  orbits  in  some  neighborhood  of  an  isolated  singular  point  the  following 
concept  is  useful: 


Definition:  Let 

^  =  A£  +  G(£)  (IV  .2.9) 

be  an  ordinary  differential  equation  as  in  (IV.2.5).  An  orbit  £  =  £(£)  is 
quasiradial  if: 

lim  £(f)  =  0.  (IV.2.10) 

t— »;fcoo 

Equivalently  we  can  express  this  definition  as: 

lim  (x(t),y(t))  =  0  lim  {p{t),q{t))  =  0  (IV.2.11) 

t-»  +  oo  t— »-foo 

lim  (z(i),  y(t))  =  0  lim  (p(l),  9(0)  =  (IV.2.12) 

t  — » —  oo  1  — * —  oo 


In  other  words,  if  the  characteristic  curves  are  quasiradial,  they  are  quasiradial  in  the 
image  plane  and  in  gradient  space.  Some  quasiradial  orbits  are  depicted  in  figure  11. 

The  stable  manifold  theorem  [ABMA80,  p.527]  and  [HAR64,  p.242]  defines  condi¬ 
tions  under  which  the  solutions  to  the  linearized  characteristic  equations  (IV. 2. 8)  have 


where  A  is  the  four  by  four  matrix: 

(0  0  2  0\ 

0  0  0  2  1 
2a  0  0  Oj 

0  27  0  0 ) 

and  where  G  has  the  following  properites: 

i)  G(t)  e  c 2 

2)  G(0)  =  0 

3)  f(0)  =  0 

Every  solution  £  =  £(t)  to  (IV.2.5)  is  called  an  orbit.  The  equation: 

is  called  the  linearized  characteristic  equation. 
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Figure  11.  Quasiradial  orbits 


the  same  topological  structure  as  the  solutions  to  the  characteristic  equations  (IV. 2.5). 
One  of  the  crucial  constraints  is  that  all  eigenvalues  of  A  have  non-vanishing  real  part. 
A  precise  definition  of  a  manifold  can  be  found  in  [ST069,  p.203|.  Here  it  suffices  to 
observe  that  a  surface  defined  by  a  Ck  function  z  =  z(x,  y)  is  a  two-dimensional  Ck 
manifold.  Before  discussing  the  stable  manifold  theorem  we  state  a  theorem  concern¬ 
ing  two  ordinary  differential  equations  in  two  unknowns  in  some  neighborhood  of  an 
isolated  critical  point: 

Theorem:  (Node  Theorem  [HAR64,  p.213j)  Let 

^  =  anx  + a12y  +  fi[x,y)  (IV.2.13) 

^  =d21x  +  a22y  +  h[x,y ) 

be  a  system  of  two  ordinary  differential  equations  where  a,j  for  i,j  =  1,2  are 
constants  such  that  an<i22  —  012021  7^  0  and  fi  for  i  =  1,  2  are  real  continuous 
functions  for  which  the  following  conditions  hold  in  some  neighborhood  of  the 
origin: 

h  =  o(|a;|  +  M)  fi  =  o{ \x\  -f  |y|)  (IV.2.14) 

Thus  (1  ,y)  =  (0,0)  is  an  isolated  critical  point.  If  both  eigenvalues  of  the 
linearized  equation  of  (IV.2.13)  are  real  and  have  the  same  sign,  then  all  orbits 
are  quasiradial  and  the  critical  point  is  called  a  node.  In  particular  if  both 
eigenvalues  are  negative,  then  each  orbit  tends  to  the  origin  as  t  — ►  -foo  and  if 
both  eigenvalues  are  positive,  then  each  orbits  tends  to  the  origin  as  t  — ►  — 00. 
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The  previous  theorem  can  be  viewed  as  a  special  case  of  the  following  theorem: 
Theorem  (Stable  Manifold  Theorem):  Let 

ft=At  +  G(  0  (IV.2.15) 

be  an  ordinary  differential  equation  where  A  is  a  constant  matrix  which  has 
di  eigenvalues  with  negative  real  part  and  d !2  eigenvalues  with  positive  real 
part,  where  di  and  d 2  are  both  positive,  and  where  G(£)  satisfies  the  following 
conditions: 

1)  <?(£)€  G2 

2)  G(0)  =  0  (IV.2.16) 

3)  If  (0)  =  0. 

Then  there  exist  an  e  >  0  such  that  (IV.2.15)  has  solutions  £  =  £(£)  ^  0 
satisfying: 

ll£(0llc<t  0  as  t  -+  00  (rV.2.17) 

and  has  solutions  £  =  £(£)  7^  0  satisfying: 

||£(t)||«“  *  ►  0  as  t  —  —00  (IV.2.18) 

Furthermore,  for  sufficiently  small  e  >  0,  the  point  £  =  0  and  the  set  of  points 
£(£)  which  satisfy  (IV. 2. 17)  sweep  out  a  unique  C 2  manifold  of  dimension  dl, 
the  stable  manifold.  Similarly,  for  sufficiently  small  t  >  0,  the  point  £  =  0  and 
the  set  of  points  £(£)  which  satisfy  (IV.2.18)  sweep  out  a  unique  C2  manifold 
of  dimension  d2,  the  unstable  manifold. 

Thus  the  curves  which  sweep  out  the  stable  (unstable)  manifold  are  quasiradial. 

IV. 3.  Integral  Surfaces  near  a  Singular  Point 

We  will  now  prove  the  main  results  of  this  chapter. 

Lemma:  Let  p2  +  q2  =  E(x,y )  be  a  constrained  eikonal  equation.  If  a  locally 
convex  solution  exists  in  some  neighborhood  of  the  singular  point,  then  it  is 
swept  out  by  quasiradial  characteristic  curves. 


Proof:  Suppose  a  locally  convex  solution  z  =  z(x,y )  exists.  As  z  =  z(x,  y)  is  assumed 
to  be  C2,  we  can  write  p  and  q  in  some  neighborhood  of  the  singular  point  as: 

p  —  aux  +  a12y  +  o(|x|  -f  |y|) 
q  =  a12x  +  a22y  4-  o(|x|  +  |y|) 


(IV.3.1) 


where  atJ  for  i,j  =  1,2  are  constants.  As  the  origin  is  a  singular  point,  p  and  q  have 
no  constant  terms.  Note  that  the  Gaussian  curvature  K  of  z  —  z(x,  y)  at  the  origin  is 
defined  by: 

K  —  011022  —  (IV.3.2) 

Recall  (section  III. 2. 4)  that  a  surface  is  locally  convex  at  a  point  if  its  Gaussian  curvature 
is  positive  there.  Substituting  the  expressions  (IV.3.1)  into  the  first  two  characteristic 
equations  (IV.2.3)  of  a  constrained  eikonal  equation  gives: 

dx 

—  =  2(anx  -f  a12y)  +  o(|x|  +  |y|)  (IV.3.3) 

~  =  2(a12x  +  a22y)  4-  o(\x\  4-  |yj). 

Using  the  node  theorem  we  deduce  that  the  characteristic  curves  in  the  x-y  plane  are 
quasiradial  if  and  only  if  both  eigenvalues  of  the  linearized  equations  have  the  same 
sign.  A  simple  calculation  shows  that  this  is  the  case  only  when  K  >  0.  Similarly,  we 
can  show  that  the  characteristic  curves  in  gradient  space  are  quasiradial  if  and  only  if 
K  >  0  which  in  turn  is  true  only  if  an  and  a22  have  the  same  sign.  Assuming  that 
K  >  0,  the  sign  of  the  eigenvalues  is  the  same  as  the  sign  of  an.  | 

In  the  next  section  we  will  actually  compute  the  coefficients  a,.,  for  i,j  =  1,2  such 
that  K  >  0. 

Finally  we  are  ready  to  state  and  prove  the  main  theorem  of  this  chapter: 

Theorem:  Let  p2  +  q2  =  E(x,y)  be  a  constrained  eikonal  equation.  Then  there 
exists  a  unique  locally  convex  solution  in  some  neighborhood  of  the  singular 
point. 


Proof:  It  follows  from  the  previous  lemma  that  if  a  locally  convex  solution  to  a  con¬ 
strained  eikonal  equation  exists,  it  is  swept  out  by  quasiradial  characteristic  curves.  So 
we  have  to  show  that  such  a  solution  exists  and  is  unique.  This  is  achieved  by  showing 
that  the  unstable  manifold  is  the  locally  convex  solution.  To  this  end  we  investigate  the 
linearized  characteristic  equations  (IV.2.8)  of  a  constrained  eikonal  equation.  An  easy 
calculation  shows  that  the  matrix  A  has  two  positive,  real  eigenvalues  and  two  negative, 
real  eigenvalues.  Thus  we  can  apply  the  stable  manifold  theorem,  which  states  that 
there  exist  exactly  two  C2  manifolds  which  are  swept  out  by  quasiradial  characteristic 
curves.  A  locally  convex  solution  therefore  exists.  From  the  node  theorem  we  get  that 
the  solution  z  —  z(x,  y)  satisfying  the  positivity  condition  is  the  unstable  manifold, 
whereas  the  stable  manifold  is  the  surface  defined  by  z  —  — z(x,y).  Hence  the  locally 
convex  solution  is  unique.  | 
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IV. 4.  Formal  Power  Series  Solution 


In  the  previous  section  we  proved  that  there  exists  a  unique,  locally  convex  solution 
to  a  constrained  eikonal  equation  in  some  neighborhood  of  the  singular  point,  but  did 
not  show  how  to  actually  compute  such  a  solution.  We  will  do  so  in  this  section  for  the 
case  where  E(x,y)  is  analytic.  In  particular,  we  shall  construct  a  formal  power  series 
solution.  In  the  case  where  the  eikonal  equation  is  not  analytic,  such  a  solution  tells 
us  only  about  the  behavior  of  an  integral  surface  of  the  equation. 


Definition:  A  formal  power  series  is  an  expression  of  the  form: 


/-*o+5>  (IV  .4.1) 

where  20  denotes  a  constant  and  for  each  j,  z3  are  homogeneous  polynomials 
of  degree  j.  We  write  the  first  order  partial  derivatives  of  /  as: 


df 

dx 


-Eft 

3  —  1 


df 

dy 


OO 


-E* 


(IV.  4.2) 


where  for  each  j,  pj  and  qj  are  homogeneous  polynomials  of  degree  j.  We  will 
say  that  p2  4"?2  —  E[x,y)  is  a  super  constrained  eikonal  equation  if  E(x,  y)  is 
C°° .  Then  /  is  a  formal  power  series  solution  to  a  super  constrained  eikonal 
equation  if: 


2 


r«/,2 


Cff  +  Cff  =  az2  4-  fry  +  IV2  + 

y  j=i 


(IV.4.3) 


where  for  each  j,  e:  is  a  homogeneous  polynomial  of  degree  j. 


Suppose  for  a  moment  that  we  have  computed  a  formal  power  series  solution  to  a 
constrained  eikonal  equation.  A  theorem  due  to  Borel  [HAR,  p.261]  states  that  there 
exists  a  C°°  function  z[x,y)  which  has  a  given  power  series  as  its  Taylor  series  expansion. 
Unfortunately  /  does  not  determine  z{x,y )  uniquely  since  two  different  functions  can 
have  the  same  power  series  expansion.  For  example  the  functions  g[x,  y)  and  g(x,y) 
defined  by  (IV. 4.4)  have  the  same  power  series  expansion  at  the  origin: 

OO  OO 

ff(*»y)  =  21 x*  +  12  y3 

3  =  1  3  =  1 

g{x>y)  —  9ix<y)  + 


(IV.4.4) 
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Note  that  the  first  order  partial  derivatives  of  g{x,y )  and  g(x,y )  are  the  same  only  at 
the  origin.  Thus  they  cannot  both  be  solutions  to  a  given  eikonal  equation.  As  a  formal 
power  series  solution  satisfies  the  eikonal  equation  at  the  origin  but  not  necessarily  at 
any  other  point,  it  is  not  necessarily  a  solution  to  a  given  constrained  eikonal  equation. 
Yet  the  formal  power  series  solution  does  tell  us,  qualitatively,  about  the  behavior  of  a 
solution  to  an  eikonal  equation.  Suppose  a  solution  z  to  the  eikonal  equation  exists  in 
some  neighborhood  of  the  origin.  Then  z  and  z  are  tangent  to  at  least  second  order  at 
the  origin. 

Lemma:  Let  p2  -\-q2  =  E(x,  y)  be  a  constrained  eikonal  equation  where  E(x,  y) 

is  analytic.  Then  its  formal  power  series  solution  is  the  solution  to  the  equation. 


Proof:  A  version  of  the  stable  manifold  theorem  proves  that  if  E(x,  y)  is  analytic,  then 
the  stable  (unstable)  manifold  is  analytic  [COLE55,  p.330j.  The  lemma  follows.  | 

Lemma:  Let  p2  -f*  q2  =  E(x,y)  be  a  super  constrained  eikonal  equation. 
Then  there  exists  a  unique,  locally  convex  formal  power  series  solution  to  this 
equation  in  some  neighborhood  of  the  singular  point. 


Proof:  (outline)  Equating  the  appropriate  terms  in  (IV. 4.3)  we  obtain 

•  an  equation  for  the  quadratic  terms  and 

•  a  recurrence  relation  for  each  of  the  higher  order  terms. 

First  (section  IV. 4.1)  we  shall  prove  that  there  is  a  unique  solution  to  the  equation  for 
the  quadratic  terms  if  we  impose  the  constraints  that  the  formal  power  scries  solution 
be  positive  and  convex.  The  next  step  (section  IV. 4. 2)  is  to  determine  the  higher  order 
terms  which  is  done  by  inductively  solving  the  recurrence  relation.  If  the  quadratic 
terms  have  been  determined  such  that  the  formal  power  series  solution  is  convex,  each 
step  of  this  induction  can  be  carried  out  uniquely.  The  first  of  the  following  two 
equations  determines  the  quadratic  terms  and  the  second  is  the  recurrence  relation 
from  which  the  higher  order  terms  can  be  calculated: 

Pi  +  ?i  =  ax 2  +  PXV  +  IV2  (IV.4.5) 

2piPfc  +  2<?! qk  =  gk+\  for  k  >  1  (IV.4.6) 

where  for  each  k,  gk  denotes  a  homogeneous  polynomial  of  degree  k.  Each  gk  is  easily 
computed  using  the  power  series  expansion  of  E{x,y)  as  we  will  show  in  section  IV. 4. 2. 

I 


IV.4.1.  Quadratic  Terms 

In  this  section  we  determine  pi  and  qi  for  given  values  of  a,  ft  and  7  in  (IV.4.5). 
Note  that  these  terms  define  the  Gaussian  curvature  of  a  solution  to  a  constrained 
eikonal  equation.  Since  p  1  and  q%  are  linear  homogeneous  polynomials  they  can  be 
written  as: 

Pi  =  anx  -f  a12y  (IV.4.7) 

Ql  —  A21z  -+■  ^221/ 


where  the  coefficients  a,j  for  i,j  =  1,2  are  constants.  As  p^  and  qy  satisfy  equation 
(IV.4.5),  the  coefficients  a.ij  for  i,j  =  1,2  are  constrained  by: 

[aux  +  a12y)2  +  [a2ix  -f  a22y)2  =  ax2  +  (3xy  -f  7y2.  (IV.4.8) 

Furthermore  px  and  q\  are  the  linear  terms  of  the  first  order  partial  derivatives  of  a 
smooth  function.  Thus  the  so-called  compatibility  condition  has  to  hold: 

dp  _  dq 

dy  dx 

which  constrains  the  coefficients  a12  and  021: 

®12  =021*  (IV.4.10) 

Equating  appropriate  terms  in  (IV.4.8)  and  using  the  compatibility  condition  we  derive 
three  equations  for  an,ai2  and  022: 

2  1  2 

all  +«12  =a 
2ai2(an  -|-  a22)  =  P 

a12  +  ®22  =  7- 

This  system  of  equations  only  has  a  solution  if  both  a  and  7  are  not  negative,  which 
was  assumed. 

As  we  wish  to  show  that  there  exists  only  one  positive  convex  formal  power  series 
solution  z(x,  y)  to  an  eikonal  equation,  we  want  the  coefficients  au  for  i,j  =  1,2  to 
satisfy  the  previous  equations  and  K  to  be  positive: 

K  —  diia22  —  ai2-  (IV. 4. 12) 


(IV.4.11) 


(IV.4.9) 


Thus  the  Gaussian  curvature  of  z  at  the  origin  is  determined  only  by  p  1  and  <71.  We  will 
say  that  a  solution  for  the  coefficients  atJ  for  i,j  —  1,2  is  convex  when  K  is  positive. 
Note  that  for  K  to  be  positive  it  is  necessary  that  a n  and  022  have  the  same  sign. 
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Figure  12.  Vectors  px  and  9, 


We  can  view  p\  and  qi  as  vectors  in  the  x-y  plane.  Then  the  vector  product  of  the 
vectors  px  and  91  is: 

Pi  X  9i  =  ana22  —  a*2.  (IV.4.13) 

Recall  that  the  vector  product  of  two  vectors  can  be  written  also  as: 

Pi  X  <71  =  |Pi||<7i|  sin  p  (IV. 4. 14) 

where  <p  denotes  the  angle  between  the  vectors  p\  and  91  as  shown  in  figure  12.  Thus 
for  pi  and  qi  to  define  a  convex  surface,  the  angle  between  them  has  to  be  less  than 
180  degrees.  We  show  now  that  the  equations  (IV. 4. 11)  constrain  this  angle.  The  inner 
product  of  pi  and  qi  is: 

Pi  •  9i  =  012(011  -f  022)-  (IV. 4. 15) 

Equivalently,  we  can  write  the  inner  product  of  pi  and  91  as: 


Pi  •  9i  =  |Pi  Iki  |  cos  ^  (IV.4.16) 

where  ip  again  denotes  the  angle  between  the  vectors  pi  and  9*.  The  following  equations 
also  hold: 


|pi|  =  l\/fl?i+0?al  (IV-417) 

|<?l!  =  I  \/  a12  +  °22l" 

So,  we  express  the  cosine  of  the  angle  p  between  p  1  and  q  1  as: 


0 

cos  P  =  — ; — ; - - 

2|>/«7| 


(IV. 4. 18) 


and  there  are  two  different  angles  y?i  and  <P2  =  27r  —  <p i  which  satisfy  this  equation. 
The  next  equation  defines  sin  <p  in  terms  of  al2,a  and  7: 


sin  <p  = 


(±V/Q~  QnXzbxA  ~  °ia)  —  0*2 

Iv^Tl 


(IV.4.19) 


Inasmuch  as  we  are  interested  solely  in  convex  solutions  (i.e.,  <p  <  180°),  this  equation 
has  to  be  solved  only  for  the  case  where  the  product  of  the  square  roots  is  positive,  or 
equivalently,  when  an  and  022  have  the  same  sign.  Combining  the  last  two  equations 
yields: 

s/a-a\2S/^~~a\2-a\2  f  /F~ 


'a  anVl  <*12  fl12  _  L  P2 

\y/cn\  V  4a^ 


(IV.4.20) 


Since  the  the  coefficients  a,  0  and  7  define  a  positive  form,  the  previous  equation  is 
well  defined  and  can  be  rewritten  as: 


4(\/4a7  —  /?2  +  a  +  7)a22  =  /?2 


(IV.4.21) 


It  is  possible  to  determine  ai2  in  terms  of  a,0  and  7  from  this  equation  as  long  as 
0^7  and  0  0.  The  case  where  0  =  7  and  /?  =  0  will  be  investigated  later. 

Case  l)a^7  and  0  0 

Let  us  denote  the  two  possible  solutions  for  012  obtained  from  (IV. 2. 21)  by  an  and 
— a12.  Without  loss  of  generality  we  assume  that  0  >  0.  The  results  for  0  <  0  are 
analogous.  The  two  solutions  for  an, 012  and  022  are: 


ail  =  V®  —  S12  ®12  =  Ol2  022  —yl— -Oj2 


(IV.4.22) 


an  =  —\jct  —  a\2  ax2  =  —an  a22  =  —\Jl  —  *n-  (IV.4.23) 

The  coefficients  defined  by  (IV.4.22)  determine  the  unique  positive  convex  solution. 
Case  2)  a  =  7  and  0  =  0 

In  this  case  equations  (IV.4.11)  can  be  written  as: 


aii  +o?2  =  « 
2ai2(an  022)  =  0 
0 12  +  0^2  =  a- 


(IV.4.24) 

(IV.4.25) 

(IV.4.26) 
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Subtracting  (IV.4.26)  from  (IV.4.24)  we  get: 

an  —  Oj2  =  0  (IY.4.27) 

and  deduce  that: 

flu  —  ifl22-  (rV.4.28) 

When  an  =  —0221  K  is  negative.  So  we  only  have  to  investigate  the  case  where 
an  =  fl22-  It  follows  from  (IV.4.25)  that  a^  =  0,  thus  we  can  express  the  coefficients 
a.ij  for  t,  j  —  1,2  in  terms  of  a  as: 


on  =  ±\/a 

ai2  =  0  (IV.4.29) 

a22  =  ~4 ~\/q. 

Both  solutions  for  the  coefficients  atJ  for  i,j  =  1,2  are  convex,  but  only  one  of  them 
is  positive.  Thus  we  have  shown  that  in  the  case  where  a  =  7  and  0  —  0  a  unique 
positive  convex  solution  exists.  | 


Actually  we  can  also  show  the  following  theorem: 

Theorem:  Let  p2-f  g2  =  E(x,  y)  be  a  constrained  eikonal  equation.  Then  at  the 
singular  point,  the  Gaussian  curvature  of  each  integral  surface  has  the  same 
absolute  value  and  is  determined  by  the  (limited)  Taylor  series  expansion  of 
E( x,  y )  at  that  singular  point. 


Proof:  Recall  that  the  curvature  at  the  origin,  denoted  by  K ,  is: 

K  =  flua22  —  °?2-  (IV. 4. 30) 

Using  equations  (IV.4.11)  we  derive  an  expression  for  this  curvature  in  terms  of  a,0,i 
and  012: 

(on  -f-  a,2i)2  — 

4°12 

(an  4*  G22)2  =  a  —  aj2  +  7  —  a?2  “1“  ^ana22  (IV.4.31) 

(an  +  a22)2  =  a  +  7  -f  2K. 

The  desired  expression  for  K  is: 

K  =  \\-fr  -(“  +  1)1-  (1V.4.32) 

l  1ai2 
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From  equations  (IV.4.18)  and  (IV.4.19)  we  get  the  two  solutions  for  a\2  in  terms  of  a,  /? 
and  7: 


(IV.4.33) 


Substituting  these  two  expressions  for  a12  into  equation  (IV.  4.32)  gives: 


(IV.4.34) 


Thus  the  absolute  value  of  the  curvature  at  the  singular  point  of  all  integral  surfaces  of 
a  constrained  eikonal  equation  can  be  directly  determined  from  the  image  intensities 
defined  by  E(x,y).  | 


IV.4.2.  Higher  Order  Terms 

In  this  section  we  calculate  the  higher  order  terms  in  a  formal  power  series  solution 
to  a  constrained  eikonal  equation  using  the  solutions  for  p\  and  qit  and  show  that  this 
can  be  done  in  a  unique  fashion  if  pi  and  qi  determine  a  locally  convex  solution.  First 
we  derive  the  recurrence  relation  which  has  to  be  solved  to  determine  the  higher  order 
terms.  Suppose  that  px  and  qi  are  known.  For  P2  and  q2  the  following  equation  must 
hold: 

2pipa  +  2qi92  =  h-  (IV.4.35) 

For  P3  and  q$  the  following  equation  holds: 

2piP3  +  2<jig3  +pl  +  ql  =  U-  (IV.4.36) 

Now,  by  assumption  p2  and  q2  are  homogeneous  polynomials  of  degree  two.  We  can 
therefore  write  the  previous  equation  as: 

2PiP3  “I"  2qi<73  =  <74  (IV.4.37) 

where  g4  can  be  determined  from  U,p2  and  q2.  Thus  to  determine,  for  each  k,  Pk  and 
qk  we  have  to  solve  the  following  recurrence  relation: 

PiPfc-i  +  9i9fc— 1  —  9k  (IV. 4. 38) 

where  gk  is  a  known  homogeneous  polynomial  of  degree  k. 

We  now  introduce  a  new  coordinate  system  denoted  by  X  and  Y : 


X  =  Px  Y  —  qi. 


(IV. 4.39) 
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We  can  do  so,  as  pi  and  qi  are  linearly  independent  vectors.  (An  easy  calculation  shows 
that  if  pi  and  qi  are  linearly  dependent  they  define  a  surface  which  has  zero  curvature 
at  the  origin.  In  other  words  the  coefficients  a,  (3  and  7  would  not  define  a  positive 
form.)  Thus  for  each  fc  >  2,  we  can  write  (IV.4.38)  as: 

XP{X,Y)  +  YQ{X,Y)  =  F{X,Y)  (IV.4.40) 

where  F(X,Y)  is  a  homogeneous  polynomial  of  degree  k  and  P(X,Y)  and  Q(X,Y)  are 
homogeneous  polynomials  of  degree  k  —  1.  This  follows  immediately  from  the  definition 
of  X  and  Y .  This  last  equation  always  has  a  solution,  i.e.,: 

P(X,Y)  =  T-"--  (IV  .4.41) 

Q{X,y)=Wl 

since,  using  these  two  functions  in  (IV.4.40),  we  obtain: 

XFx{X,  Y)  +  YFy{X,  Y)  =  kF{X,  Y)  (IV.4.42) 

which  is  Euler’s  homogenity  relation  [BRSE75,  p.246]  for  a  homogeneous  polynomial 
of  degree  k. 

Now  we  show  that  any  solution  for  P{X,Y)  and  Q(X,Y )  can  be  written  as: 

P(X,  Y)  =  P{X,  Y)  -f  C(X,  Y)Y  (IV.4.43) 

Q(X,  Y)  =  Q(X,  Y)  -  C(X,  Y)X 

where  C(X,Y)  is  a  homogeneous  polynomial  of  degree  k  —  2.  Suppose  that  Pl{X,Y), 
Q1(X,Y )  and  P2(X,Y),Q2(X,Y )  satisfy  (IV.4.40).  Then  P(X,Y)  =  P1  —  P2  and 
Q(X,  Y)  —  Q1  —  Q 2  have  to  satisfy  the  following  equation: 

xp(x,y)  +  yQ(x,y)==o.  (iv.4.44) 

We  are  looking  for  two  functions  P(X,Y)  and  Q(X,Y)  which  satisfy  the  last  equation 
and  are  homogeneous  polynomials  of  degree  k  —  1.  Furthermore  P(X,YJ  must  vanish 
at  Y  =  0  whereas  Q(X,Y)  must  vanish  at  X  =  0.  Thus  P(X,Y)  and  Q(X ,Y)  must 
be  of  the  form: 


P{X,Y)=  C(X,Y)Y  (IV.4.45) 

Q(X,Y)=~C(X,Y)X 


which  proves  our  assertion  (IV.4.43). 
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Hence  we  can  write  the  functions  P(X,Y)  and  Q(X,Y),  which  satisfy  equation 
(IV.4.40),  as: 


P(X,  Y)  =  +  C{X,  Y)Y  (IV.4.46) 

K> 

Q(X,Y)  =  -  C{X,Y)X. 


Yet  we  want  to  find  some  P(X,Y)  and  Q(X,Y)  which  not  only  satisfy  (IV.4.40)  but 
also  the  compatibility  condition: 


dP_dQ 
dy  dx 


(IV.4.47) 


We  will  show  that  if  pi  and  q\  define  a  convex  solution,  only  one  homogeneous  poly¬ 
nomial  C(X,Y )  of  degree  k  —  2  exists  such  that  P  and  Q  satisfy  this  condition.  Using 
(IV.4.46),  the  partial  derivatives  of  P  and  Q  with  respect  to  x  and  y  are: 


%  -  fr1 + + r(§f + if> 

f = j  <F-  S + F-  f  >  -  f 1 f + i  £>• 


(IV.4.48) 


Recall  the  definitions  of  X  and  Y  in  terras  of  x  and  y: 


X  =  anx  -f-  ai2y  (IV.4.49) 

Y  =  al2x  +  a22y. 


So  the  partial  derivatives  of  X  and  Y  with  respect  to  x  and  y  are: 


ax 


dY 

dx 


ai2 


ax 

dy 

dY 

dy 


ai2 


a22- 


(IV.4.50) 


As  F{X,Y)  is  a  homogeneous  polynomial,  Fxy  —  Pyx ■  Using  equations  (IV.4.48)  and 
(IV.4.50)  we  can  express  the  compatibility  condition  (IV.4.47)  as: 


•^[PxxO’U  —  FyY^j.2  +  Pxy(o 22  ~  Oil)]  =  (IV.4.51) 

—  C(X,Y)(a„  +  022)  —  Cx(Xau  +  Y  012)  —  Cy{Xc  12  +  Y  022)- 


We  wish  to  determine  when  there  is  only  one  polynomial  C(X,Y)  of  degree  k  —  2 
which  satisfies  the  previous  equation.  In  so  doing,  we  will  compute  the  coefficients  in 
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the  polynomial  C(X  ,Y).  Then  let 

k 

F(X,Y)  =  £'f>.k-,X'Yk-' 

«= o 

k 

fx(x,y)  = 

i=l 

Fr(X,  Y)  =  jrj  /<ifc_,(fc  -  t)X‘rfc— 1 
1=0 

Fxx{X,Y )  =  'Srf,,k-.ii(i-l)X‘-'Yk-‘ 

i=2 

Fyy(X,Y )  =  -  0(*  -  *  -  l)X*rfc-‘~2 

1=0 

/■xr(x.r)  =  f  ~‘)x-‘rk—\ 

1=1 

Similarly,  letting 

<?(x,r)  -  k’E^.k-i-^Y‘-’-k 

i=o 

c*(x,y)  =  flci.k-i-tixi-'Y"-*-2 

7  =  1 

Cr(X,y)  =  £cilfc_,-_2(fc  -  j  -  2)x>rfc-'-3 

7—0 

gives: 

xcx(x,y)  =  f;ciffc-i_2y^rfc-'-2 

7=1 

k—2 

YCx(X,Y)=  J2c),k-j- ijXi-'Y*-*-1 

7  =  1 

XCr(X,y)  =  J2c:,k-:-2{k  -  j  -  2)X’  +  1Yk ~’~3 

7=0 

k—3 

YCy(X.Y)  =  X)  -  j  -  2)X‘Yk-‘-‘. 

7=0 


(IV.4.52) 


(IV.4.53) 


(IV.4.54) 
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Using  the  last  three  systems  of  equations  in  (IV.4.51)  we  can  rewrite  the  compatibility 
condition  as: 


-  l)Xi-2rfc-‘a12  - 

*  i=2 
k—2 

£  -  *X*  -  *  -  i)x‘r*— 2ai2  4 

t=0 

k— 1 

'52fi,k-ii(k  -  -  au)]  = 

i=i 

k—2 

-  53  cJ,fc_J_2XJVfc— ■J— 2(an  +  a22)  - 
i=o 

fc— 2 

53ci,fc_i_2y^Vfc-J'-2o11- 

k—2 

i=l 
fc— 3 

53  c,-  fc_y-2(*  -  j  -  2)X>+lYk-’-*al2  - 

y=o 

fc— 3 

53  cif  *_,■_,(*  -  j  -  2 )X’  Yk-*-aa2a. 
y= o 


(IV.4.55) 


Equating  the  coefficients  of  the  powers  XttYk  **  for  /x  =  0, A:  —  2  we  obtain  the 
following  equations: 


For  n  =  0: 

^(/2,fc— 22ai2  —  fo,kk(k  —  l)a  12  4  fi,k—i{k  —  l)(a22  —  an)]  =  (IV.4.56) 

—  co,fc_ 2(an  4  a22)  —  ci,*— 3012  —  c0ik-2{k  —  2)a22 

For  0  <  /i  <  ft  —  2: 

^[/M4-2,fc-M-2(M  4-  2)(/i  +  l)ai2  —  /M,fc-„(fc  —  p)(fc  —  M  —  l)ai2  4 

/,*+i,fc-M-i(fc  —  M  —  1)(*  —  l*)(a22  —  an)]  =  (IV.4.57) 

—  c/1)(c_M_2(an  4  a22)  —  cMifc_^-_2poii  — 

—  4  l)ai2  —  c/i_itfc_M_1(A:  —  p  —  l)ai2  — 

2(*  —  M  —  2)a22 
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-^[fk.oKk  —  l)an  —  fk— 2,2^12  4-  fk—  l,l(fe  —  1)(«22  ~  Oil)]  =  (IV.4.58) 

—  c*-2,o(«u  +  ^22)  —  Ck— 2,o(*  —  2)an  —  Ck— 3,1^12- 


Denote  by  FM,  for  /x  =  0, k  —  2,  the  left  hand  side  of  each  appropriate  equation, 
multiplied  by  — 1.  Then  we  can  rewrite  (IV.4.56),  (FV.4.57)  and  (IV.4.58)  as: 


Fo  — Co,fc — 2(^11  +  022(fc  —  1)1  +  Ci(fc_3ai2 

F„  =cM_lffc —  n—  1  )a12  +  (IV.4.59) 

cM,fc_M^2(an(Ai  +  1)  +  ai2{k  —  n  —  1))  + 

CM+l,fc— M— 3(M  “I-  1)^12 

^*ie — 2  =CAc— 3, 1^12  +  Ck_2,o[au(fc  ~  1)  +  a22]- 


Now,  for  each  fixed  A  >  2  we  need  to  determine  the  coefficients  c, 2  from  the 
previous  system  of  &  —  1  linear  equations.  This  system  can  be  written  in  matrix  form 
as: 

Ak^ic  =  F  (IV.4.60) 

and  has  a  unique  solution  if  the  matrix  Afc_i  is  nonsingular  (i.e.,  its  determinant  is 
not  equal  to  zero).  Note  that  Ak— 1  is  a  tridiagonal  band  matrix,  so  its  determinant 
can  be  computed  using  the  recurrence  relation  (IV.4.61)  which  we  now  derive. 

Let  Bn  be  an  n  X  n  matrix,  let  Bn— 1  denote  the  submatrix  of  Bn  obtained  from 
Bn  by  eliminating  its  top  row  and  leftmost  column  and  let  Bn—2  be  similarly  obtained 
from  Bn_  1.  Then  Bn  can  be  written  as: 


fa  b 

Bn=  c  Bn_! 

VO 


(IV.4.61) 


The  determinant  of  Bn  may  be  computed  in  terms  of  the  determinants  of  its  sub- 
matrices: 

detBn  =  adetB„_i  —  bcdetBn_2.  (IV.4.62) 

It  suffices  to  prove  by  induction  that  Ak— 1  is  nonsingular  if  an  022  —  <*12  >  0. 
Without  loss  of  generality  (see  section  IV. 4.1)  we  can  assume  that  an  >  0  and  022  >  0. 
The  fith  row  of  Ak— 1  is: 


(k  —  fi  —  l)aj2  Qii(m  +  1)  “I"  <i22[k  —  fi —  1)  ai2(^  +  !)•  (IV.4.63) 

Set  j  =  k  —  fi  —  2  and  denote  by  A j  the  determinant  of  the  submatrix  of  Ak— 1 
consisting  of  its  j  1  bottom  rows  and  j  -j- 1  rightmost  columns.  We  prove  that  Ak— 1 
is  nonsingular  by  induction  on  j. 
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For  the  basis  case,  we  have  to  show  that  Ao  7^  0  and  Ai  7^  0: 

Ao  =  dn{k  —  1)  4"a22  (IV. 4. 64) 

Ai  =  [au[k  —  2)  -f"  2a22][fln(&  —  1)  +  022]  —  —  2)* 

By  assumption  Ao  is  positive.  Rewrite  Ax  as: 

Ai  =  a2n{k  -2 ){k  -  1)  +  2a11a22(A:  -  1)  +  2a*2  +  (aua22  -  a212)(k  -  2).  (IV.4.65) 

Since  we  have  assumed  that  ana22  —  a22  >  0,  it  follows  that  Aj  is  always  positive. 
Assume  now  that  Ay_2  and  Ay_x,  for  j  >  1,  are  positive.  We  show  that  A y  is 


then  also  positive  which  completes  the  proof.  Using  (IV.4.62)  we  get: 

A  y  =  [on(fc  —  j  —  1)  +  a22(i  +  l)]A,_i  —  a\2(k  —  j  —  1);  Ay_2.  (IV.4.66) 

where  Ay_i  is  defined  by: 

Ay — 1  =  [an(A:  —  j)  -f  a22>]Ay_2  —  a22(k  —  j){j  —  l)Ay _3  (IV.4.67) 

and  where  we  define  A_i  to  be  1.  Using  (IV.4.67)  in  (IV.4.66)  gives: 

Ay  =au{k  —  j  —  l){(«u(fc  —  j)  -F  o22y]Ay_2  —  a22(k  —  j)(j  —  l)Ay_3}  -f 

<222(7  4  l)Ay_x  —  a\2(k  — j  —  l)y'Ay — 2  (IV .4.68) 

Ay  =(ana22  —  ah){k  —  j  —  l)jAy_2  -j-  022(7  +  l)Ay_i  + 

«n(*  —  7  —  l)(fc  —  7)[au Ay— 2  —  0*2(7  ~  1)^j— s]*  (IV.4.69) 

To  verify  that  Ay  is  positive,  one  need  only  show,  therefore,  that: 

all Ay — 2  >  <*12(7  —  l)Ay_3.  (IV.4.70) 

We  proceed,  once  again,  by  induction.  For  7  ==  3  we  wish  to  show: 

ouAi  >  2a22A0.  (IV.4.71) 


Using  (IV.4.64)  and  (IV.4.65)  in  (IV.4.71)  this  can  be  equivalently  written  as: 

a  1 1  [a  11  —  2)(fc  —  1)  4*  2aii022(fc  —  1)  4  2022  4*  (011^22  —  ai2)(^  —  2)]  >  (IV.2.72) 

2gj2[oii(&  —  I)  4*  022]. 

To  see  that  the  previous  equation  holds  it  suffices  to  observe  that: 

2a?,a22(*  —  1)  >  2a22an{k  —  1) 

2ana22  >  20j2fl22' 


(IV.4.73) 
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This  solves  the  basis  case.  Assume,  by  induction,  that  the  next  inequality  holds: 


anAj— 2  >  a\2(j  —  1)Aj_3. 

(IV.4.74) 

We  need  to  show  that: 

anAj_1  >  a\2j£^j—2- 

(IV. 4.75) 

Using  (IV.4.67)  in  this  last  inequality  gives: 

aii(^  —  j)Aj— 2  +  aiifl22jAj_2  —  ajiai2(A:  —  j)[j  —  l)Ay_3  >  a\2j A;_2.  (IY.4.76) 

This  holds  if: 

ah(k  —  y)A>-2  >  anal2(/c  —  j)(j  —  1)Aj_3 

(IV.2.77) 

which  is  true  by  the  induction  hypothesis  (IV.4.74). 

Hence  there  exists  a  unique  locally  convex  formal  power  series 
strained  eikonal  equation  which  can  be  easily  computed.  | 

solution  to  a  con- 
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Chapter  V 


B-Silhouettes 


V.l.  Overview  of  Basic  Concepts 

As  discussed,  our  goal  is  to  find  necessary  constraints  such  that  an  image  can  be 
interpreted  in  a  unique  way  when  its  image  irradiance  equation  is  known.  It  was  shown 
in  chapter  III  that,  in  general,  there  can  be  an  infinite  number  of  different  surfaces 
which  satisfy  the  same  image  irradiance  equation;  in  other  words,  the  image  of  each 
of  these  surfaces  is  the  same  for  a  fixed  imaging  configuration.  Recall  that  in  sections 
III. 2. 2  and  III.3  we  discussed  how  to  detect  a  b-silhouette  in  an  image.  It  remains 
now  to  investigate  whether  the  existence  of  such  b-silhouettes  can  be  used  to  interpret 
an  image.  In  this  chapter  we  will  identify  three  constraints  upon  an  image  irradiance 
equation,  one  upon  the  reflectance  map,  one  upon  the  o-silhouette  and  one  upon  the 
function  E(x,y).  If  these  constraints  hold  for  some  image  irradiance  equation,  only  one 
surface  defined  by  a  C2  function  which  satisfies  the  equation  exists. 

What  kind  of  information  about  an  integral  surface  can  one  deduce  from  a  singular 
image  irradiance  equation?  Let  R[p,q)  =  E(x,y)  be  a  fixed,  singular  image  irradiance 
equation  whose  nondegenerate  b-silhouette  is  defined  by  w(x,y)  =  0.  As  previously 
discussed  (section  III. 3),  the  surface  normal  to  the  bounding  contour  of  an  integral 
surface  is  parallel  to  {±wx,  ;^u;v,  0). 


V.2.  Uniqueness  Theorem 

In  this  section  we  obtain  constraints  which  assure  that  if  an  image  irradiance 
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equation  has  a  C2  integral  surface,  it  is  unique.  So  let 

R(p,q)  =  E{x,y)  (V.2.1) 

be  a  singular  image  irradiance  equation.  Consider  the  following  constraints  upon  this 
equation: 

Cl)  R{p,q)  =  p2  +  q2. 

C2)  The  b-silhouette  defined  by  w(x,  y)  =  0  is  a  closed,  smooth  curve  in  the  x-y 
plane.  Furthermore,  the  points  ( x,y )  at  which  the  image  irradiance  equation  i3 
defined,  lie  in  the  region  bounded  by  this  b-silhouette. 

C3)  The  function  E(x,  y)  has  exactly  one  stationary  point  [xq,  y0)  and  satisfies 
the  following  conditions  in  some  neighborhood  of  (zo,  Vo):  E(x o,yo)  =  0, 
E[x,y )  >  0  for  ( x,y )  7^  (z0,yo)  and  E[x,y)  vanishes  precisely  to  second  order 
at  (z0,  yo). 


Uniqueness  Theorem:  Let  R(p,q)  =  E[x,y)  be  an  image  irradiance  equation 
for  which  constraints  Cl,  C2  and  C3  hold  and  suppose  a  C 2  integral  surface 
defined  by  z  —  z(x,  y)  of  this  equation  exists.  Then  the  only  other  solution  to 
the  equation  is  z  =  — z(x,y). 


Proof:  Let  R(p,q)  ==  E(x,y)  be  a  fixed  image  irradiance  equation  for  which  constraints 
Cl,  C2  and  C3  hold.  First  note  that  the  point  P  —  ( x,y,p,q )  =  (xo,yo,0, 0)  is  an 
isolated  singular  point  of  the  image  irradiance  equation  (see  also  section  IV. 2).  There 
are  then  two  observations  which  allow  us  to  prove  the  theorem.  First,  as  the  b-silhouette 
is  a  closed  curve,  an  integral  surface  of  the  equation  has  to  be  compact.  Second,  from 
constraints  Cl  and  C3  we  can  deduce  that  such  a  surface  is  convex  at  the  singular 
point,  which  allows  us  to  apply  results  of  the  previous  chapter. 

Suppose  z  —  z(x,y)  defines  a  C2  integral  surface  of  an  image  irradiance  equation 
as  defined  in  the  uniqueness  theorem.  Then  from  C2  we  may  infer  that  z  defines  a 
compact  surface  (section  III. 2. 4).  Note  also  that  z(x,y)  is  defined  for  every  point  (x,y) 
which  lies  within  or  on  the  b-silhouette  and  therefore  has  a  bounding  contour.  Thus 
there  exists  a  point  P  at  which  z  has  an  extremum.  In  particular  the  tangent  plane  at 
P  is  parallel  to  the  x-y  plane. 

From  condition  C3  we  can  deduce  that  there  exists  a  point  (x0,y0,zo)  on  2  such 
that  the  plane  tangent  to  z  at  this  point  is  parallel  to  the  x-y  plane,  i.e.,  p(xo,yo)  =  0 
and  q(x 0,  yo)  =  0.  By  assumption  this  is  the  only  singular  point  of  the  image  irradiance 
equation  and  therefore  the  only  point  on  z  where  the  tangent  plane  is  parallel  to  the 
x-y  plane.  Furthermore,  as  the  image  irradiance  equation  is  singular,  the  point  (zo,yo) 
lies  in  the  interior  of  the  region  of  the  x-y  plane  bounded  by  the  b-silhouette.  Hence 
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P  =  (x0,y0,zQ).  By  the  assumption  on  E{x,y),  z  is  either  convex  at  P  or  has  a  saddle 
point.  Since  P  is  the  point  where  the  surface  has  maximal  (minimal)  height,  z  must  be 
convex  there. 

In  chapter  IV  we  proved  that  if  an  image  irradiance  equation  satisfies  Cl  and  C3, 
there  exists  only  one  positive  and  exactly  one  negative  convex  solution  denoted  by  zp 
and  zn  =  — Zp  respectively.  Thus  there  are  exactly  two  integral  surfaces  z  —  z{ x,y) 
and  z  =  —  z(x,  y).  | 

By  using  transformation  methods,  we  can  enlarge  the  class  of  singular  image  irradiance 
equations  for  which  the  uniqueness  theorem  holds.  Let 

/( V  +  2  Bpq  +  C92  +  2  Dp  +  2  Eq)  =  E{x,  y )  (V.2.2) 

be  a  singular  image  irradiance  equation  where  /  is  a  bijection  and  A,  B,  C,  D  and  E 
are  real  constants  such  that  6  >  0  and  AS  <  0  where  6,  A  and  5  are  defined  in  the 
following  equations: 

6  =  AC  —  B2 
A  B  D 

A=  B  C  E  (V.2.3) 

DEO 
S  =  A  +  C. 

The  constraints  upon  the  constants  A,B,C,D  and  E  in  equation  (V.2.2)  assure  that 
the  curves  R{p,q)  —  c,  for  any  constant  c,  are  closed.  Let  the  b-silhouette  of  the 
equation  be  a  closed  and  smooth  curve.  Then  (V.2.2)  can  be  transformed  into  an  image 
irradiance  equation  of  the  form  (V.2.1)  for  which  C2,  holds  as  is  shown  in  appendix 
II.  If,  after  the  transformation,  the  function  E[x,y)  satisfies  C3,  then  the  uniqueness 
theorem  holds  for  (V.2.2). 

The  next  corollary  follows  directly  from  the  uniqueness  theorem  in  this  chapter. 
We  will  abbreviate  \J  %2  +  y2  by  r. 

Corollary:  Let  p2  -f-  q2  =  E[r)  be  an  image  irradiance  equation  where  E(r) 
satisfies  constraints  C2  and  C3.  Suppose  a  C 2  integral  surface  of  this  image  ir¬ 
radiance  equation  exists.  Then  it  is  rotationally  symmetric  and  can  be  obtained 
by  integrating  E(r).  In  this  case  the  b-silhouette  is  a  circle. 

Proof:  First  we  write  the  eikonal  equation  in  polar  coordinates: 

zr  +  ^4  =  E(r)- 


(V.2.4) 
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Let  z  =  z(r )  define  the  rotationally  symmetric  integral  surface  of  the  above  eikonai 
equation.  Thus  i,j(r)  =  0  and  we  can  compute  both  rotationally  symmetric  solutions  by 
integrating  ±\/E(r).  It  follows  from  the  uniqueness  theorem  that  the  image  irradiance 
equation  has  only  rotationally  symmetric  integral  surfaces.  | 


Note  that  the  above  corollary  does  not  hold  if  the  image  does  not  contain  a 
b-silhouette.  In  section  III. 2. 5  we  showed  that  the  integral  surfaces  of  a  continuous 
rotationally  symmetric  eikonai  equation  are  not  themselves  necessarily  rotationally 
symmetric.  The  uniqueness  result  for  the  special  image  irradiance  equation  (III.3.4) 
has  been  independently  obtained  by  [DSY80]. 


V.3.  Counterexamples 

In  the  previous  section  we  discussed  sufficient  constraints  under  which  the  solution 
to  a  singular  image  irradiance  equation  is  unique.  Are  these  constraints  necessary? 
Although  we  are  not  able  to  answer  this  question  completely,  we  now  shed  some  light 
upon  it.  In  particular  we  try  to  find  the  class  of  image  irradiance  equations  for  which 
most  likely  there  is  no  set  of  constraints  that  assure  the  existence  of  only  one  global 
solution. 

Image  irradiance  equations  satisfying  the  constraints  of  the  uniqueness  theorem 
have  closed  iso-brightness  curves,  i.e.,  the  curves  R(p,  q)  —  c  are  closed.  So  let  us 
examine  singular  image  irradiance  equations  whose  iso- brightness  curves  are  not  closed. 
Such  an  image  irradiance  equation  is  given  by: 


P  +  9  = 


— (»  +  y) 

y/l  —  (x2  +  y2) 


(V.3.1) 


While  constraint  C2  holds  for  (V.3.1),  an  image  irradiance  equation  where  the  reflectance 
map  is  R(p,q)  =  p  +  q  never  has  a  singular  point.  The  general  solution  to  (V.3.1)  is: 


z 


(x,  y)  =  ^1  —  (z2  +  y2)  +  w(y  —  x ) 


(V.3.2) 


where  w  is  any  C 1  function.  Figures  13  through  16  illustrate  some  solutions  to  equation 
(IV.3.1). 

Another  example  of  an  image  irradiance  equation  whose  iso-brightness  curves  are 
not  closed  is  given  by: 


_  xy 

Pq  l-(x2  +  y2)' 


(V.3.3) 


This  equation  satisfies  constraint  C2,  i.e.,  its  b-silhouette  is  a  closed  and  smooth  curve. 
Furthermore  (0,0)  is  the  singular  point  of  (V.3.3)  and  E(x,y )  vanishes  precisely  to 
second  order  there  although  E{x,y)  is  not  positive  in  the  neighborhood  of  (0,0).  One 
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of  the  solutions  to  (V.3.3)  is  the  sphere: 


z{x,y)  =  J 1  —  (i2  +  !/2) 


(V  .3.4) 


whereas  another  surface  satisfying  (V.3.3)  is: 


z(x,  y)  =  f(t)  +  x3  —  y2  where 
t  =  1  —  (x2  +  y2)  and 

m = t\lh+l + l[in(\/  h + 1 +i)~  1°l'/if+'  ~ 1"' 


(V.3.5) 


Recall  that  coastraiut  C2  expressed  the  fact  that  the  b-silhouette  is  a  closed  and 
smooth  curve.  We  now  demonstrate  that  if  the  b-silhouette  does  not  obey  C2  our 
uniqueness  result  does  not  hold.  Previously,  it  was  not  expected  that  the  umquenes 
results  for  image  irradiance  equations  containing  closed  b-silhouettes  would  be  differen 
from  the  Results  for  those  containing  open  b-silhouettes.  An  example  of  an  equation 
for  which  Cl  holds,  but  whose  b-silhouette  is  not  a  closed  curve  is. 


+  ?2  =  tz  +  l 


(V.3.6) 


Equation  (V.3.6)  does  not  have  a  singular  point.  A  solution  to  this  equation  is. 


z{x,y)  =  +  y 


(V.3.7) 


i 


Figure  17.  z(x,  y)  =  y/x  +  y 


which  is  shown  in  figure  17.  Two  other  solutions  to  (V.3.6)  are: 

*(*,„)  =  V  5  +  '>  +  iN  \/h  +  1  + *>  _  ln(\/s 


+  1  -  1)|.  (V.3.8) 


y/x[l  —  8z)  x  n 

z(x>  y )  = - ^ - atan  I/  a - 1  +  \/3  y 

2  a./o  V  8x 


2  4n/2 

A  more  complicated  counterexample  is: 

P  +<7  = 


a:2  +  y2 

1  —  xy 


(V.3.9) 


This  equation  clearly  satisfies  Cl  but  its  b-silhouette  is  not  a  closed  curve.  The  function 
E{x,  y)  vanishes  precisely  to  second  order  at  the  singular  point  of  (V.3.9)  and  is  positive 
in  the  neighborhood  of  the  origin.  Two  different  solutions  to  (V.3.9)  are  given  by 
equation  (V.3.10)  (shown  in  figure  18)  and  (V.3.11)  (shown  in  figure  19). 


z{x,y)  —  2\/l  — 


xy 


z{x,y)  =  (i  —  xy)\ 


1  —  xy 


—  1  —  atam 


1  —  xy 


-1  + 


x 2  —  y 2 


(V.3.10) 

(V.3.11) 


r 


f 
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j 

Figure  18.  z(x,y)  =  2y/l  —  xy  | 

Note  that  the  surface  defined  by  (V.3.11)  is  C1  at  the  origin  (although  not  C2)  and  is 
only  defined  in  two  quadrants,  whereas  the  image  irradiance  equation  is  defined  in  all 
four  quadrants.  The  surface  defined  by  equation  (V.3.10)  is  not  convex  at  the  origin. 

Only  when  the  b-silhouette  is  a  closed  curve  can  we  deduce  that  a  surface  is  convex 
at  the  singular  point,  an  observation  which  allows  us  to  prove  the  uniqueness  theorem. 

For  the  following  eikonal  equation  we  give  two  solutions  which  are  both  C2  at  the  origin 
but  only  one  of  which  is  convex: 

-j 

► 


Equation  (V.3.12)  has  a  singular  point  and  E(x,  y)  satisfies  constraint  C 3.  Two  solutions 
to  this  equation  are: 


z(x,y) 

z(*>  y) 


arcsin(xy) 

\/l  —  x2y2  -|- 


(V.3.13) 
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Chapter  VI 


Numerical  Methods 


VI. 1.  General  Remarks 

In  this  chapter  we  discuss  some  numerical  methods  proposed  to  solve  the  shape 
from  shading  problem.  The  first  shape  from  shading  algorithm  was  implemented  by 
Horn  [HO70],  who  proposed  to  solve  the  characteristic  equations  (III. 2. 2)  using  standard 
numerical  methods  for  solving  ordinary  differential  equations.  Horn  observes  that  the 
surface  gradient  at  a  singular  point  P  (section  III. 2. 2  and  chapter  IV)  is  uniquely  defined 
by  the  image  irradiance  equation.  By  assuming  that  a  surface  is  locally  convex  in  some 
neighborhood  of  P  and  then  estimating  the  curvature  of  this  surface,  he  is  able  to 
calculate  an  initial  curve  (sections  A.6  and  A.7).  Horn  notes  that  using  this  heuristic, 
only  one  surface  is  calculated.  However,  it  is  not  guaranteed  that  there  exists  a  solution 
to  any  given  image  irradiance  equation  which  is  locally  convex  in  some  neighborhood  of 
P  and  so  the  surface  determined  using  Horn’s  algorithm  may  not  be  an  integral  surface. 
Furthermore,  the  algorithm  computes  at  most  one  of  the  oo.  integral  surfaces  of 
an  image  irradiance  equation.  We  proved  in  chapter  IV  •  v  i»_  ;  case  of  an  eikonal 

equation  where  some  technical  conditions  are  imposed  on  l£[x,  y )  a  unique  locally  convex 
solution  exists  in  some  neighborhood  of  a  singular  point.  In  this  situation  then,  Horn’s 
algorithm  computes  the  unique  convex  solution.  (Note  that  using  our  results  it  is  not 
necessary  to  estimate  the  curvature  at  a  singular  point  in  order  to  uniquely  compute 
the  surface.)  Unfortunately,  the  algorithm  is  slow,  numerically  unstable,  and  relies  on 
the  presence  of  singular  points. 
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Strat  [STR79]  developed  an  iterative  shape  from  shading  algorithm  which  we 
discuss  in  detail  in  section  VI. 2.  One  of  the  major  shortcomings  of  this  algorithm 
is  that  an  initial  curve  lying  on  a  surface  is  required  as  an  input. 

Brooks  [BR079]  suggests  a  Waltz-like  [WA75]  constraint  propagation  scheme  to 
determine  the  shape  of  a  surface  from  its  shading.  Like  Strat,  Brooks  assumes  that  a 
surface  is  C2  and  that  an  initial  curve  which  is  embedded  in  the  surface  is  known. 
In  addition,  he  imposes  an  upper  limit  upon  the  curvature  of  a  surface  that  can 
be  determined.  This  constraint  (which  stems  from  experiments  with  human  visual 
systems)  is  needed  for  the  algorithm  to  converge. 

The  three  methods  mentioned  above  can  only  solve  the  shape  from  shading  problem 
under  the  assumption  that  a  surface  is  smooth  and  does  not  have  a  bounding  contour. 
If  an  image  contains  a  b-silhouette  it  can  be  analyzed  by  using  the  algorithm  due  to 
Ikeuchi  and  Horn  [IKH081]  which  is  explained  in  section  VI.3. 


VI. 2.  St  rat’s  Algorithm 

Strat  [STR79]  proposes  an  iterative  algorithm  to  solve  an  image  irradiance  equation 
under  the  assumption  that  an  integral  surface  is  defined  by  a  C2  function  z  =  z(x,y). 
He  assumes  that  image  irradiance  is  measured  at  discrete  points  on  the  image  plane 
and  that  the  reflectance  map  is  given.  In  this  scheme,  a  square  grid  is  imposed  on  the 
x-y  plane.  We  denote  a  grid  point  by  the  tuple  (i,j),  the  image  irradiance  measured  at 
(i,  j)  by  Eltj,  the  first  order  partial  derivatives  of  z  —  z(x,y )  (which  defines  an  integral 
surface)  at  a  point  ( i,j )  by  pi>3  and  qx>3,  and  the  reflectance  map  evaluated  at  [Pi,j,qx,j) 
by  Rtj .  The  objective  of  the  algorithm  is  to  calculate  pM  and  qi<3  at  every  point  (t,j). 
Strat  does  not  address  the  problem  of  how  to  compute  the  function  which  defines  the 
surface,  i.e.,  zXi3,  from  the  values  of  pXi3  and  qxj. 

So  let  R{p,q)  =  E(x,y)  be  a  given  image  irradiance  equation.  We  now  show  how 
to  obtain  a  function  e  (called  the  error  function )  in  the  variables  pt>J  and  qx<3  which  is 
zero  when  the  image  irradiance  equation  is  satisfied  at  every  point  and  pXt3  and  qii3  are 
the  first  order  partial  derivatives  of  z  =  z( x,  y).  Strat’s  algorithm  computes  iteratively 
the  values  for  px<3  and  qx<3  which  minimize  the  function  e.  This  function  is  the  weighted 
sum  (over  all  (t,  j))  of  two  functions  ej  •  and  e"3  whose  derivations  we  now  explain.  The 
function  ej  •,  defined  by  (VI. 2.1),  is  zero  when  p,i3  and  qxj  satisfy  the  image  irradiance 
equation: 

=  (Btj  -  Hi,,)2.  (VT.2.1) 

The  values  for  p,i3  and  qXiJ  which  minimize  et  cannot,  however,  be  chosen  inde¬ 
pendently;  they  are  the  first  order  partial  derivatives  of  a  C2  function  z  =  z{x,y). 
Consequently,  equation  (VI. 2. 2)  holds  for  p  and  q  [DIR72]: 


Py  =  Qx- 


(VI. 2. 2) 
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GJ  +  i) 


OJ) 


(i  +  1j+ t) 


(i  +  1j) 


Figure  20.  Template 


Strat  assumes  that  the  image  irradiance  equation  is  solved  over  a  simply  connected 
region  in  the  x-y  plane.  Hence  (VI. 2.2)  can  also  be  expressed  in  integral  form  [DIR72]: 

j>  pdx  -j-  qdy  =  0.  (VI.2.3) 

We  now  show  how  a  discrete  approximation  of  equation  (VI.2.3)  is  used  to  derive  e'y. 
Strat  suggests  the  loop  around  four  adjacent  points  on  the  grid  as  a  path  over  which 
the  integral  in  (VI.2.3)  should  be  taken.  Such  a  path  is  depicted  in  figure  20.  A  discrete, 
first  order  approximation  to  equation  (VI.2.3)  is  then  given  by: 

P*,3  “1“  Pi+i,j  ~f~  9»-t-i ,j  4"  9«+i,j+ 1  P»+i,y-fi  Pt,j  + 1  9t,y-f  l  Qij  —  0.  (VI. 2. 4) 

Thus  the  values  for  ptij  and  qtiJ  which  the  algorithm  calculates  must  not  only  satisfy 
the  image  irradiance  equation,  but  also  equation  (VI. 2. 4).  In  general,  a  point  (i,  j) 
belongs  to  four  templates.  For  each  template,  a  discrete  approximation  of  (VI.2.3)  can 
be  derived  and  ptij  and  qij  must  satisfy  each  such  approximation.  So  e*  •  is  defined 
as: 

€i,j  —{Pi+ij  4-  q<+i,j  4-  qi+ i.y+1  —  p<+i.i+i  — 

Pi, 3  + 1  —  9i,i+i  +  Pi,  j  —  9t,i)2  + 

(Pm— i  4  Pi+I,j~ l  +  i,i— i  4  9«+i,i  ~ 

P*  +  1,J  —  9«,j  — 1  —  Pi,j  —  qi,j )2  + 

(Pt-l.j  — 1  -f  Pi,j—  1  4-  qi,j  —  1  —  Pi—  1,3  — 
qt-i,j  —  ?t— i,i— i  —  Pi, 3  4-  g.j)2  4- 

{Pi—i,j  4  Qi,3  + 1  Pi,i  +  1  Pi  —  1,J  + 1 

—  9t— i,j  4-  Pi,j  4  9«,j)2- 


(VI. 2. 5) 
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The  total  (global)  error  e  is  then  defined  as: 

e  =  £«,+K>)  (vn.e) 


where  p  is  a  constant  scalar  chosen  appropriately  to  equate  the  dimensions  and  to 
weigh  the  magnitudes  of  the  errors  and  ■.  At  each  iteration,  Strat’s  algorithm 
minimizes  e.  Hence  at  every  point  (i,  j)  the  following  two  equations  must  be  solved  for 
Pi,:  and  qi<3: 


de 

dPiJ 

de 

dqi.j 


=  0 


=  0. 


(VI. 2. 7) 


The  Gauss  Seidel  method  is  used  to  solve  equations  (VI. 2.7). 

It  remains  to  describe  how  the  initial  assignment  for  the  values  p,^  and  qli3  is 
chosen.  Strat  assumes  that  the  algorithm  is  applied  to  solve  an  image  irradiance 
equation  over  a  rectangular  region  in  the  x-y  plane.  For  points  (r,  j)  on  the  boundary 
of  the  rectangle,  the  values  for  p,tJ  and  qhj  are  given  as  input  data  and  for  interior 
points,  pij  and  qt|J  are  set  to  zero.  As  mentioned  in  section  VI.  1,  Strat’s  method  is  not 
very  useful  in  practice,  since  it  requires,  as  input,  the  initial  values  for  ptiJ  and  at 
the  boundary  of  the  domain  in  which  the  algorithm  is  applied.  It  is  not  shown  whether 
the  solution  computed  by  the  algorithm  is  independent  of  this  initial  choice  of  values 
for  pitJ  and  qtiJ  at  the  interior  points. 

Experimental  results  appear  to  support  the  conjecture  that  Strat’s  algorithm  con¬ 
verges,  although  a  proof  of  convergence  is  not  given.  One  of  the  shortcomings  of  Strat’s 
algorithm  is  that  its  performance  depends  upon  the  order  in  which  the  grid  points  are 
scanned,  i.e.,  the  order  in  which  the  values  for  pit:  and  q,tj  are  updated. 


VI. 3.  Ikeuchi  and  Horn’s  Algorithm 

Ikeuchi  and  Horn  [IKH081]  designed  and  implemented  an  iterative  algorithm  which 
differs  in  various  respects  from  the  ones  described  in  the  previous  sections. 

•  Their  assumption  about  the  surface  whose  shading  is  to  be  analyzed  is  weaker. 
They  observe  that  for  a  surface  to  look  smooth  it  suffices  for  the  function 
defining  it  to  be  only  C1.  Recall  that  under  this  assumption,  an  integral  surface 
of  an  image  irradiance  equation  can  be  built  from  characteristic  strips  (see 
chapter  III  and  appendix  I). 

•  They  show  that  a  b-silhouette  can  be  used  as  initial  data,  which  is  not  possible 
using  Brooks’  or  Strat’s  algorithm. 
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One  of  the  strongest  features  of  the  method  is  that  it  incorporates  data  obtained  from 
a  b-silhouette.  Recall  that  the  surface  normal  at  a  point  on  the  bounding  contour  is 
parallel  to  the  normal  vector  to  the  b-silhouette.  Before  proceeding  to  describe  and 
analyze  Ikeuchi  and  Horn’s  algorithm,  we  review  their  slightly  different  formulation  of 
an  image  irradiance  equation. 

In  chapter  III  it  was  shown  that  an  image  irradiance  equation  defines  a  relation 
between  the  image  plane  and  gradient  space.  In  the  original  formulation  of  this  equation 
[H075]  the  ( p,q )  coordinate  system  was  chosen  to  specify  gradient  space.  With  this 
underlying  coordinate  system  an  image  irradiance  equation  can  be  formulated  as  a 
FOPDE.  On  the  other  hand  this  system  has  the  disadvantage  that  for  points  on  the 
bounding  contour  p  and/or  q  assume  infinite  value.  In  order  to  find  a  numerical  solution 
to  an  image  irradiance  equation  it  would  be  useful  to  transform  gradient  space  into 
a  different  space  SP,  where  (u,  v)  is  the  underlying  coordinate  system,  so  that  the 
components  of  a  surface  normal  vector  (which  always  have  finite  value),  instead  of  p 
and  q,  can  be  computed.  To  achieve  this  we  have  to  find  a  space  SP  and  a  mapping 
m  between  gradient  space  and  SP,  satisfying  the  following  constraints: 

•  The  mapping  m  between  gradient  space  and  a  region  of  SP  is  one-to-one  and 
onto. 

•  The  mapping  m  maps  the  whole  (unbounded)  p-q  plane  into  a  bounded  region 
in  SP. 

Recall  (section  III. 2. 5)  that  gradient  space  can  be  obtained  by  projection  of  the  Gaussian 
hemisphere  from  its  center  onto  a  plane  placed  at  its  south  pole.  Ikeuchi  and  Horn 
observe  that  by  choosing  a  different  projection  of  the  Gaussian  hemisphere  onto  a  plane, 
each  point  on  this  hemisphere  gets  mapped  into  a  point  other  than  infinity  on  the  plane. 
Several  projections  are  suggested,  one  of  them  being  the  stereographic  projection.  We 
can  think  of  this  projection  in  geometric  terms  as  a  projection  onto  a  plane  tangent  to 
the  hemisphere  at  the  south  pole.  The  center  of  projection  is  the  north  pole,  not  the 
center  of  the  sphere.  In  this  case  the  resulting  coordinate  system,  denoted  by  ( f,g ),  is 
defined  in  terms  of  p  and  q  by: 

f  =  +p2  +  g2  — 1) 

p 2  -f  q2 

2q{y/1  +  Pi  -t-g2  ~  1) 

9  p2  -f  q2 

We  can  now  write  an  image  irradiance  equation  in  the  form: 

R{f,g)  =  E(x,y)-  (VI-3.3) 

Next  we  describe  Ikeuchi  and  Horn’s  method  of  solving  the  shape  from  shading 
problem.  The  basic  concepts  of  this  algorithm  are  similar  to  those  described  in  sec¬ 
tion  VI. 2.  The  inputs  to  this  algorithm  are  the  measured  image  irradiance  and  the 


(VT.3.1) 

(VI.3.2) 
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reflectance  map  R(f,g)  and  the  desired  outputs  are  the  values  fXtJ  and  gx>j.  (The 
problem  of  how  to  determine  Zji7  from  the  values  of  /t  J  and  gxo  is  not  addressed  by 
Ikeuchi  and  Horn.)  There  are  two  major  differences  between  this  and  Strat’s  procedure. 
First,  a  different  error  function  is  used,  and  second  the  initial  values  for  /f|J  and  r,t7 
are  computed  for  points  on  the  b-silhouette.  (The  case  where  there  are  no  b-silhouettes 
in  the  image  is  discussed  later.)  For  points  not  on  the  b-silhouette,  /t>J  and  gtiJ  are 
initialized  to  zero  just  as  in  Strat’s  algorithm. 

We  now  explain  Ikeuchi  and  Horn’s  derivation  of  the  error  function.  They  exploit 
the  assumption  that  the  function  z  =  z[ x,y)  defining  an  integral  surface  is  C1.  A 
function  F  =  F{x,y)  is  continuous  at  a  point  (zo,2/o)  if,  for  any  e  >  0,  there  exists  a  6 
such  that,  for  any  (i,  y)  which  lies  in  the  circle  defined  by  (i  —  Xq)2  -)-  (y  —  y0)2  <  6, 


we  have: 

\F{x,  y)  —  F(x0,  y0)|  <  (VI.3.4) 

A  discrete  approximation  of  equation  (VI.3.4)  is  given  by: 

l-Fi-fi.j  +  P'i — i,y  -j-  F%,j+ 1  +  Fi,j— i  —  4F<,y|  <  4e.  (VI.3.5) 

The  error  function  e  is  therefore  defined  to  be: 

e  —  ^(e* ,  +  XeV)  where  (VI.3.6) 

~  Ri,if  (VI-3-7) 

ei,j  +  fi,j- 1-1  +  —  4 fi,j)2  + 

(ff«-4-i,y  +  4*  9i,j+i  4-  9t,j— l  —  4g,,y)2-  (VI.3.8) 


In  equation  (VI.3.6)  X  denotes  a  scalar  factor  chosen  appropriately  to  equate  the 
dimensions  and  to  weigh  the  magnitudes  of  the  errors  ej  •  and  ej  •. 

As  in  Strat’s  algorithm,  the  function  e  is  minimized.  Thus,  at  every  iteration  the 
next  two  equations  are  solved  for  f,j  and  y,i7: 

=  0  (VI.3.9) 

°u,j 


^9\,j 

Once  again  the  Gauss  Seidel  method  is  used  to  solve  equations  (VI.3.9).  A  lower  bound 
on  the  number  of  iterations  needed  using  this  method  is  proportional  to  the  mesh  size 
but  an  upper  bound  is  not  known.  However,  numerical  experiments  indicate  that  indeed 
the  algorithm  converges. 

As  mentioned  before,  f,t]  and  gtiJ  are  initialized  to  zero  at  every  point  except  those 
on  the  b-silhoii'’,te.  It  is  not  shown  that  the  algorithm  is  independent  of  these  initial 
values. 

In  the  case  where  there  are  no  b-silhouettes  in  the  image,  heuristics  have  to  be 
used  to  initialize  the  algorithm.  Ikeuchi  and  Horn  do  not  specify  the  number  of  points 
at  which  the  exact  values  for  fltJ  and  gltJ  need  to  be  known  in  order  to  guarantee  that 
the  algorithm  will  compute  an  answer. 
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Chapter  VII 


Conclusion 


In  this  report,  we  have  investigated  the  question  of  how  much  information  con¬ 
cerning  the  shape  of  an  object  can  be  deduced  from  its  shaded  image.  Even  assuming 
that  adequate  data  is  available  to  derive  an  image  irradiance  equation  is  insufficient  to 
solve  the  reconstruction  problem  uniquely;  in  general,  for  a  fixed  imaging  configuration 
there  are  many  surfaces  which  have  the  same  shaded  image.  Thus  our  goal  has  been 
to  identify  constraints  by  which  the  reconstruction  problem  can  be  solved  uniquely. 

We  first  analyzed  the  continuous  image  irradiance  equation.  The  information 
needed  to  restrict  the  solutions  to  such  an  equation  to  a  single  one,  was  previously 
known  [COHI62b]:  If  a  strip  is  specified  which  is  not  a  characteristic  strip  and  along 
which  A  0  (III. 2. 3),  then  there  exists  a  unique  surface  which  contains  this  strip 
and  whose  image  has  a  particular  shading  for  a  fixed  imaging  configuration.  We  were 
interested  in  whether  edges  could  constitute  an  initial  curve.  In  the  case  where  the 
image  irradiance  equation  is  linear,  it  is  not  possible  to  distinguish  among  the  multiple 
surfaces  which  could  give  rise  to  a  known  edge.  For  nonlinear  image  irradiance  equations 
an  edge  constrains  the  possible  surfaces  to  a  small  number.  However,  if  there  exists  a 
surface  which  contains  a  vertex  and  edges  emanating  from  it,  it  is  unique. 

Then  we  discussed  how  singular  points  of  an  eikonal  equation  constrain  its  possible 
solutions.  In  particular  we  proved  that  there  exists  a  unique  (up  to  translation  in 
the  a-direction)  positive  convex  surface  which  satisfies  an  eikonal  equation  in  some 
neighborhood  of  a  singular  point. 

Finally  we  investigated  images  in  which  b-silhouettes  can  be  detected.  An  image 
irradiance  equation  which  describes  the  relationship  between  the  gradient  of  a  surface 
whose  image  contains  a  b-silhouette,  and  the  shading  of  this  surface,  can  be  singular. 
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We  showed  that  singular  image  irradiance  equations  may  be  solved  using  the  same 
methods  applied  to  find  the  integral  surfaces  of  a  continuous  image  irradiance  equation. 
However,  our  ultimate  ambition  was  to  answer  the  following  question: 

Is  there  a  set  of  constraints  which  assure  that  if  an  image  irradiance  equation 
has  a  solution,  it  is  unique? 

We  answered  this  question  affirmatively  in  chapter  V.  It  was  shown  there  that  if  three 
constraints  are  known  to  hold,  the  information  about  the  imaging  situation  and  the 
surface  as  captured  by  an  image  irradiance  equation,  allow  one  to  reconstruct  the  shape 
of  the  surface  in  a  unique  manner.  Furthermore  one  can  easily  check  whether  an  image 
irradiance  equation  satisfies  these  constraints.  It  is  surprising  that  our  uniqueness 
theorem  holds  only  when  the  b-silhouette  is  a  closed  curve  (constraint  C2). 

In  order  to  evaluate  the  usefulness  of  our  uniqueness  theorem,  we  need  to  know 
which  of  the  commonly  arising  image  irradiance  equations  actually  obey  the  above 
mentioned  restrictions.  In  his  paper  on  hill-shading  [II079],  Horn  discusses  eighteen 
different  reflectance  maps  which  are  applied  to  solve  that  problem.  Constraint  Cl  holds 
for  five  of  those  reflectance  maps.  These  equations  are  of  the  form: 

R{p,q)  =  f{p2 +  q2)  (VII.1.1) 

where  /  is  a  bijection.  A  reflectance  map  of  the  form  (VII.1.1)  describes,  for  instance, 
the  situation  when  the  object  is  Lambertian  and  the  light  source  and  the  viewer  have 
the  same  position.  Furthermore  eikonal  equations  can  be  also  used  to  automatically 
analyze  images  taken  by  a  scanning  electron  microscope  which  seems  to  be  one  of  the 
prime  applications  of  our  uniqueness  theorem. 

There  are  several  issues  not  treated  here  which  would  be  interesting  to  investigate 
further:  mutual  illumination,  shadows  and  specularities,  for  example.  Another  open 
problem  is  to  determine  the  class  of  functions  R  =  R(p,  q)  which  are  reflectance  maps. 
Our  only  restriction  upon  these  functions  has  been  that  they  be  Cx .  If  a  smaller  set 
can  be  found  which  is  known  to  contain  all  reflectance  maps,  better  methods  could 
be  developed  to  solve  the  shape  from  shading  problem  in  practice.  A  further  research 
problem  is  to  find  a  good  numerical  shape  from  shading  algorithm. 
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Appendix  i 


Mathematical  Details 


This  appendix  contains  a  summary  of  results  known  about  first  order  partial 
differential  equations.  A  thorough  treatment  of  this  material  can  be  found  in  [COHI62b, 
pp. 62-153]  and  [SMI68b,  pp.259-301]. 

A.1.  Basics 

For  simplicity  of  exposition,  we  will  discuss  only  first  order  partial  differential 
equations  involving  a  function  z  of  two  variables,  z  and  y.  The  results  can  be  generalized 
to  functions  of  n  variables  in  a  straightforward  fashion.  Let  p  and  q  denote  the  first 
order  partial  derivatives  of  z  with  respect  to  z  and  y  respectively.  Then  the  relation: 


F(x,y,z,p,q)  =  0 


(A.  1.1) 


where  F  is  a  function  of  z,  y,  z,  p  and  q,  is  called  a  first  order  partial  differential  equation 
(abbreviated  in  the  following  by  FOPDE).  The  relation  (A.1.1)  is  a  linear  FOPDE  if  it 
is  linear  in  p  and  q  with  coefficients  depending  only  on  z  and  y  and  (A.1.1)  is  quasi-iinear 
if  it  is  linear  in  p  and  q  with  coefficients  depending  on  z,y  and  z. 

A  function  z( z,y)  is  called  a  solution  to  (A.1.1)  if  in  some  region  of  the  z-y  plane 
the  function  and  its  derivatives  identically  satisfy  the  equation  in  z  and  y.  Such  a 
function  is  also  called  an  integral  surface. 

The  general  solution  to  a  FOPDE  is  a  whole  set  of  solutions,  each  of  which  satisfies 
the  equation.  Given  a  FOPDE,  what  kind  of  constraints  can  be  imposed  so  that  there 
is  only  one  integral  surface  which  satisGes  both  the  constraints  and  the  equation?  Such 
constraints  arc  for  example  boundary  conditions  or  initial  values.  In  section  A. 4  we  will 
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sta,te  precisely  what  comprises  a  general  solution  and  in  sections  A.6  and  A.7  discuss 
what  kind  of  constraints  are  necessary  to  pin  down  a  particular  solution. 

Unless  otherwise  stated,  •  we  assume  that  F,  z  and  all  relevant  derivatives  are 
continuous. 

A. 2.  The  Quasi-Linear  FOPDE 

We  consider  quasi-linear  FOPDE’s  first  as  their  geometric  interpretation  is  rather 
clear  and  so  the  relevant  method  for  solving  them  can  be  explained  and  understood 
easily.  In  this  case  we  rewrite  the  relation  (A.1.1)  as: 


y.  Ap  +  b(x>  y>  Ay  =  c(x>  y,  A- 

Furthermore  we  assume  that: 


(A.2.1) 

(A.2.2) 


a2  4-  b2  0. 

Suppose  that  the  solutions  to  (A.2.1)  are  written  implicitly  as: 

G(x,y,z)  =  0.  (A.2.3) 

Differentiating  (denoted  by  subscripts)  (A.2.3)  with  respect  to  x  and  y  yields: 

Gx  -f-  Gzzx  =  0  and 

or  equivalently: 


Gx  . 

p  =  — —  and 


Gy  “f*  GzZy  —  0 
Gu 


'v 

Gm 


(A.2.4) 

(A.2.5) 


Using  these  equations  in  (A.2.1)  we  obtain: 

a(x,  y,  z)Gx(x,  y,  z)  -f-  b{ x,  y,  z)Gy(x,  y,  z )  -f  c(x,  y,  z)Gz(x,  y,  z)  —  0. 


(A.2.6) 


Note  that,  in  general,  (A.2.1)  is  a  nonlinear  FOPDE  for  the  function  z{x,y)  whereas 
(A.2.6)  is  a  linear  FOPDE  for  G(x,y,  z).  We  can  interpret  the  coefficients  a,  6  and  c 
in  (A.2.6)  as  the  components  of  a  vector  field  which  is  defined  by  £  =  £{x,y,  z)  — 
\a(x,y,  z),  b(x,y,z),  c(x,y,  a)).  Then  we  can  rewrite  (A.2.6)  as: 


tf,VG)  =  0 


(A.2.7) 


where  VG  denotes  the  gradient  of  G  and  (, )  the  inner  product  of  two  vectors.  Equation 
(A.2.7)  expresses  the  fact  that  £  is  perpendicular  to  VG.  Since  the  vector  VG  is 
perpendicular  to  the  surface  defined  by  G{x,y,  z)  =  0  at  each  point  (x,  y,  2),  we  deduce 
that  £  lies  in  the  tangent  plane  of  this  integral  surface  at  that  point. 
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A  field  line  is  a  curve  whose  tangent  at  every  point  has  the  same  direction  as  the 
field  vector  there.  An  integral  surface  of  (A.2.7)  can  be  built  up  from  field  lines  (called 
characteristic  curves  in  this  context)  of  the  vector  field  £.  To  reiterate  the  previous 
statements:  The  tangent  at  each  point  of  a  characteristic  curve  has  the  same  direction 
there  as  the  vector  £  and  therefore,  by  virtue  of  (A.2.7),  is  perpendicular  to  the  surface 
normal  to  the  integral  surface  G(x,  y,  z )  =  0.  This  does  not  mean  that  each  quasi-linear 
FOPDE  has  a  single  solution.  Such  a  FOPDE  only  constrains  the  possible  orientations 
of  the  tangent  planes  at  each  point  to  a  one-parameter  manifold.  As  (A.2.1)  is  linear 
in  p  and  q  at  every  point  of  any  integral  surface,  all  feasible  tangent  planes  intersect  in 
a  line  which  is  called  the  Monge  axis. 

We  describe  now  a  method  for  finding  characteristic  curves.  Such  curves  can  be 
defined  as  functions  of  one  parameter  s:  x  =  x(s),  y  =  y(s)  and  z  =  a(s).  The  vector 
[x(s),  t/(s),  2(s)]  is  denoted  by  Then  ^  =  [^T^i  has  the  same  direction 

as  £  and  so  the  outer  product  of  and  £  must  equal  zero: 


,dz  dy 

b - c—  =  0 

ds  da 

dx  dz 

c*-“£  =  0 

=  0. 

ds  as 


The  relations  (A.2.8)  are  normally  written  as: 

dx:dy:dz  =  a:b:c. 


(A.2.8) 


(A.2.9) 


The  solutions  to  the  equations  (A.2.9)  comprise  a  two-parameter  family  of  curves 
in  space  (a  family  of  characteristic  curves),  yet  only  a  one-parameter  subset  of  them 
generates  the  solutions  to  the  FOPDE.  To  find  this  subset,  an  arbitrary  function 
between  the  two  free  parameters  obtained  when  solving  (A.2.9)  is  introduced.  Hence  the 
general  solution  to  (A.2.1)  contains  this  arbitrary  function.  So  each  surface  produced 
by  a  one-parameter  family  of  characteristic  curves  is  an  integral  surface.  Conversely 
we  now  show  that  each  integral  surface  is  generated  in  such  a  fashion. 

On  each  integral  surface  z  =  z(x,  y)  the  equations: 

^  =  a(x,  y,z)  ~~  b(x,  y,  z)  (A.2.10) 

define  a  one-parameter  family  of  curves:  x  —  x(s),y  =  y(s),z  —  z(x(s),y(s)).  Note 
that  on  any  curve  in  this  family: 

~  =  c(x,  y,  z)  (A. 2.11) 
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as:. 


dz  dz  dx  dz  dy 

Ts  =  aiTs  +  aiT,  =  ap  +  ^  ~ c ■ 


(A.2.12) 


Thus  every  integral  surface  is  swept  out  by  a  one- parameter  family  of  characteristic 
curves.  The  following  example  illustrates  these  ideas. 


Example: 

The  following  FOPDE  is  to  be  solved: 

F{x,y,z,p,q)  =  xp  +  yq  —  z  =  0.  (A.2.13) 

As  a(x,y,z)  =  x,  b[x,y,z )  =  y  and  c(x,  y,z)  =  1,  the  equations  for  the  characteristic 
curves  (A.2.9)  are: 

dx:dy:dz  =  x:y:l  (A.2.14) 

and  have  as  their  solution  the  two- parameter  curves  in  space: 

y  =  C\x  (A.2.15) 

z  =  C2x 


where  Cj  and  C2  are  constants.  Any  integral  surface  can  be  built  from  the  curves 
described  by  equations  (A.2.15)  and  every  such  surface  is  a  one-parameter  manifold  of 
these  curves.  Each  of  these  manifolds  is  determined  through  coupling  Cx  and  C2  by 
an  arbitrary  relation  w: 

-  =  C2  =  w{Cx)  =  «;(-).  (A.2.16) 

z  z 

Thus  the  solution  to  the  FOPDE  is: 

z  —  w(-)x.  (A.2.17) 

x 

Writing  this  in  parametric  form  gives: 

y  =  C\x  and  z  =  w(Ci)x.  (A.2.18) 


A. 3.  The  General  FOPDE 


"•r  i  iTrti  mAh  r  ■- 


We  can  apply  methods  similar  to  those  developed  in  the  previous  section  to  solve 
a  general  FOPDE: 

F{x,y,z,p,q)  =  0  (A. 3.1) 
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where  we  require  that: 

Fl  +  F2q^  0.  (A.3.2) 

Our  goal  is  to  transform  the  problem  of  finding  a  solution  to  (A.3.1)  into  the  problem 
of  integrating  a  set  of  ordinary  differential  equations.  Again,  geometrical  reasoning  will 
help  to  find  these  equations. 

Let  P  be  a  given  point  on  an  integral  surface.  Then  the  quantities  p  and  q,  which 
determine  the  direction  of  a  tangent  plane  at  P,  are  constrained  by  (A.3.1)  to  a  one- 
parameter  family  of  curves.  (In  other  words:  once  the  coordinates  x,  y  and  z  of  a  point 
P  are  fixed,  (A.3.1)  is  an  equation  for  p  and  q.  To  write  this  equation  in  parametric 
form  only  one  parameter  is  needed.)  The  envelope  of  the  tangent  planes  is  a  conical 
surface  which  can  have  several  sheets  and  is  called  the  Monge  cone.  (A  conical  surface 
is  produced  by  moving  a  straight  line  which  is  fixed  at  one  point,  along  a  curve.)  “The 
considerations  here  refer  merely  to  a  suitable  small  range  of  tangent  planes,  e.g.,  a 
portion  of  a  sheet  of  the  cone  where  q  can  be  expressed  as  a  single- valued  differentiable 
function  of  p”  [COHI62b,  p.75].  Each  generator  of  this  cone  represents  a  possible 
direction  of  the  tangent  plane  at  P  and  is  called  a  characteristic  direction.  Thus  the 
integral  surface  has  to  fit  into  the  field  of  Monge  cones,  i.e.,  has  to  always  be  tangent 
to  them. 

Recall  now  that  in  the  quasi-linear  case,  a  Monge  cone  degenerates  to  a  Monge 
axis.  To  determine  the  solution  in  that  case,  we  find  the  characteristic  curves  which 
at  every  point  have  as  their  tangent  direction  the  direction  of  the  Monge  axis  there. 
An  integral  surface  is  swept  out  by  these  curves.  In  the  case  of  a  genera!  FOPDE,  the 
same  basic  idea  works.  Again,  we  have  to  find  the  curves  which  at  every  point  have  as 
their  tangent  direction  a  characteristic  direction.  Let  such  a  curve  (called  focal  curve) 
be  given  by  x  =  x[s),y  —  y(s)  and  z  =  z(s).  Remember  that  a  one-parameter  family 
of  such  curves  should  sweep  out  an  integral  surface  z[x,y)  of  a  given  FOPDE;  in  other 
words,  the  functions  x(s),  y(s),  z(s),p(s)  and  q(s)  have  to  satisfy  this  FOPDE.  The  focal 
curves  determine  x(s),y(s)  and  z(s).  Yet,  (A.3.1)  gives  only  one  relationship  between  p 
and  q  and  so  in  order  to  determine  p  and  <7,  another  equation  is  needed.  This  equation 
can  be  obtained  by  requiring  that  focal  curves  be  embedded  in  an  integral  surface.  (A 
focal  curve  is  embedded  if  in  some  neighborhood  of  the  projection  of  this  curve  onto 
the  x-y  plane,  z  is  a  single-valued,  twice  continuous  differentiable  function  of  x  and 
y.)  Focal  curves  which  satisfy  this  last  condition  are  called  characteristic  curves  and  a 
one-parameter  family  of  characteristic  curves  sweeps  out  an  integral  surface. 

The  problem  which  actually  has  to  be  solved  is  that  of  finding  an  integral  surface. 
So  the  solution  to  this  proceeds  in  the  opposite  direction  from  that  described  in  the 
preceding  paragraphs.  First,  a  set  of  equations,  called  the  characteristic  equations  has 
to  be  found.  The  characteristic  curves  are  a  subset  of  the  solutions  of  this  set,  from 
which  an  integral  surface  can  be  built  up.  In  the  following  paragraphs  the  technical 
prerequisites  needed  to  find  an  integral  surface  of  (A.3.1)  are  developed. 

As  a  first  task  we  determine  the  equation  of  the  Monge  cone.  Again,  let  P  be  a 
fixed  point  which  has  the  coordinates  (x,y,z).  Then  p  and  q,  which  satisfy  (A.3.1), 
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are, written  as  functions  of  a  parameter  u  and  all  feasible  tangent  planes  at  (x,y,  z)  are 
expressed  as: 

(Z  -  z)  =  (X  -  x)p(u)  +  (Y  -  y)q(u).  (A.3.3) 

The  envelope  of  the  planes  described  by  the  previous  equation  defines  the  Monge  cone, 
a  conical  surface  whose  vertex  is  (x,  y,  z).  The  equation  of  the  Monge  cone  is  found  by 
eliminating  u  from  (A.3.3)  and  from  (A. 3. 4)  which  is  obtained  by  differentiating  (A.3.3) 
with  respect  to  u: 

0  =  (A.3.4) 

du  du  ' 

Differentiating  (A.3.1)  with  respect  to  the  parameter  u  gives: 


n  _  F  dP  ,  P 
°  ~  f’Tu  + f'7u 


(A.3.5) 


Assuming  that  neither  and  jf*  nor  Fp  and  Fq  vanish  simultaneously,  we  derive  the 
following  equation  from  (A.3.3),  (A.3.4)  and  (A.3.5): 


X  —  x  _  Y  —  y  _  Z-z 
Fp  Fq  ~~  pFp  +  qFq ' 


(A.3.6) 


By  substituting  all  possible  values  for  p  and  g  (i.e.,  all  values  for  p  and  q  which 
satisfy  (A.3.1))  we  obtain  all  generators  of  the  Monge  cone  at  the  point  (x,  y,  z).  “Space 
curves  having  a  characteristic  direction  at  each  point  shall  be  called  focal  curves” 
[COHl62b,  p.76].  Therefore  focal  curves  have  to  satisfy  the  following  differential 
equations: 

dx  _  dv  dz 

ds=Fp  •  Ts=Fq  ~fo=pFp  +  qFq- 

To  show  that  every  integral  surface  z  can  be  generated  in  such  a  fashion,  let  z  —  z(x,  y) 
be  an  integral  surface  on  which  p  and  q  is  also  known.  Then  the  equations: 


dx 

ds 


—  Fp 


(A.3.8) 


define  a  one-parameter  family  of  curves.  On  these  curves: 


dz 

ds 


dx  dy 
-p7s+qT, 


holds  and  using  (A.3.8)  in  (A. 3. 9)  we  obtain  the  third  equation  of  (A. 3.7): 


(A.3.9) 


Ta=pfr+ if,- 


(A. 3.10) 
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The  above  equation  is  called  the  strip  condition.  “It  states  that  the  functions  x(s),  y(s), 
2(s),p(s)  and  q{s)  not  only  define  a  space  curve,  but  simultaneously  a  plane  tangent  to 
it  at  every  point.  A  configuration  consisting  of  a  curve  and  a  family  of  tangent  planes 
to  this  curve  is  called  a  strip”  [COHI62b,  p.77] . 

Equation  (A.3.9)  states  that  the  curves  defined  by  (A.3.7)  are  focal  curves.  Now  it 
is  also  required  that  a  focal  curve  be  embedded  in  an  integral  surface.  Differentiating 
(A.3.1)  with  respect  to  x  and  y  we  obtain: 

FpPx  +  FqQx  +  FzP  +  Fx  =  0  (A.3.1 1) 

FpPy  +  Fqqy  -j-  Faq  -f  Fy  —  0. 


Using  (A.3.7)  and  the  fact  that  pv  =  qx  leads  to  the  following  equations: 


is  =  +  =  V'‘Fr  +  q'F’' 

Finally,  using  (A.3.11)  in  (A.3.12)  yields: 


fa+F,v  +  Fx  =  0 
—  0- 


(A.3.12) 


(A.3.13) 


In  summary:  If  a  focal  curve  is  embedded  in  an  integral  surface  then  along  this 
curve  the  coordinates  x,  y,  z  and  the  quantities  p  and  q  satisfy  the  following  five  ordinary 
differential  equations: 


dx 

ds 

dp 

ds 


FP 


~(pF,-\-Fx) 


dz 

ds 

dq 

ds 


=  pFP-\-qFq 


=  -(qF,  +  Fv). 


(A.3.14) 


We  can  consider  this  system  of  differential  equations  (which  are  called  characteristic 
equations)  by  itself,  i.e.,  disregarding  that  we  obtained  it  with  a  given  integral  surface 
in  mind.  Note  first  that  F{x,y,z,p,q)  is  constant  along  each  solution  to  the  system 
(A.3.14)  since: 


dF 

ds 


dP  4.  F  dq  4-  F  dZ  4-  F  dX  4-  F  Ay  - 

~F”Z  +  F-Z+F‘Z  +  F**  +  “ 


=  -  FpfpF,  +  f.)  -  F,{qF,  +  Fy)  +  F,(pFr  +  qF,)  +  F,Fr  +  F,F,  (A.3.15) 
=0. 
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Thus  F(x,y,z,p,q)  =  c,  where  c  is  a  constant,  is  a  solution  to  (A. 3.14).  The  sys¬ 
tem  of  characteristic  equations  defines  a  four-parameter  family  of  curves.  By  impos¬ 
ing  the  additional  constraint  that  the  solutions  to  (A.3.14)  also  satisfy  the  FOPOE 
F(x,y,  z,p,q)  =  0  we  obtain  a  three- parameter  subset  of  the  solutions.  “Every  solu¬ 
tion  of  the  characteristic  differential  equations  which  also  satisfies  the  equation  F  =  0 
will  be  called  a  characteristic  strip ;  a  space  curve  x(s),  y(s),  z(s)  bearing  such  a  strip  is 
called  a  characteristic  curve”  [COHI62b,  p.79j.  The  fact  that  a  one-parameter  subset 
of  the  three-parameter  family  of  curves  sweeps  out  the  integral  surface  has  already 
been  established.  Thus  the  problem  of  finding  a  solution  to  a  FOPDE  is  equivalent  to 
integrating  the  system  of  five  (or  equivalently  four  if  the  equations  are  not  written  in 
parametric  form)  ordinary  differential  equations  (A.3.14).  Note  that  since  the  charac¬ 
teristic  curves  depend  on  the  solution,  their  range  of  influence  cannot  be  determined  in 
advance. 

In  the  next  section  we  discuss  the  notion  of  a  complete  integral  and  then  show  how 
to  choose  the  appropriate  one- parameter  subset  of  the  solutions  to  the  characteristic 
equations. 

A. 4.  General  Solution,  Complete  and  Singular  Integral 

We  showed  in  the  previous  section  that  each  solution  to  a  general  FOPDE  is  swept 
out  by  a  one-parameter  family  of  curves.  Thus  we  can  write  the  equation  of  an  integral 
surface  as  a  function  of  the  coordinates  x  and  y  and  an  arbitrary  function  of  one 
variable.  Such  an  equation  is  called  the  general  solution  to  a  FOPDE. 

Suppose  now  for  a  moment  that  a  solution  z  =  <P(z,  y,  a,  b)  to  a  FOPDE  is  known, 
where  <P  depends  on  the  two  parameters  a  and  b.  Then  <P(x,  y,  a,  b )  is  called  a  complete 
integral  if  A,  defined  by: 

A  =  *Pxa$yb  ^xb&ya  (A. 4.1) 

is  not  equal  to  zero.  This  condition  assures  that  <P  really  depends  on  a  and  6,  i.e.,  that 
there  is  no  7  =  g{a,b)  such  that  &(x,y,a,b )  =  $(x,y,y). 

From  the  two-parameter  family  of  surfaces  defined  by  &(x,  y,  a ,  b),  we  can  chose 
a  one-parameter  subset  by  introducing  an  arbitrary  function  w  which  relates  a  and 
6,  e.g.,  by  setting  b  =  w(a).  Note  that  the  family  4>{x,  y,  a,  w(a))  is  a  solution  to  the 
FOPDE.  The  envelope  of  the  family  4>{x,y,a,w[a))  is  again  a  solution  to  the  FOPDE 
since  at  each  point  it  touches  a  member  of  the  family  &(x,  y,  a,  w(a)),  i.e.,  a  solution. 
We  obtain  the  equation  of  this  envelope  by  eliminating  the  parameter  a  from  the  two 
equations: 

z  =  <P{x,y,a,w(a))  (A.4.2) 

$a(z,  V,  a,  w{a))  +  4>b{ x,  y,  a,  to(a))u/(a)  =  0. 

We  assume  throughout  that  all  eliminations  are  possible  and  that  during  the  course 
of  this  process  only  functions  with  continuous  derivatives  are  obtained.  Eliminating  a 
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from  (A.4.1)  yields  an  expression  involving  an  arbitrary  function  w  of  one  variable 
which  is  a  solution  to  the  FOPDE  and  therefore  the  general  solution.  We  demonstrate 
this  fact  analytically  by  differentiating  the  first  equation  of  (A.4.2)  with  respect  to  x 
and  y: 

zx  =  #x  +  (&a  4-  <PbV>\a))ax  (A.4.3) 

zv  =  +  (#<*  +  $bV>'{a))ay. 

Recall  that  0(x,y,a,w(a))  is  a  solution  to  the  FOPDE  for  any  choice  of  the  parameter 
a.  Using  (A.4.2)  (i.e.,  4>a  -j-  0bU)'{a)  =  0)  ‘n  previous  equations  establishes  the  fact 
that  for  all  x  and  y  the  values  of  z,  zx  and  zy  are  the  same  as  those  of  0,0X  and  0y 

So  if  a  complete  integral  of  a  given  FOPDE  is  known,  we  can  obtain  the  general 
solution  by  differentiation  and  by  elimination  of  parameters.  (This  latter  process  can 
in  practice  be  tedious  or  impossible  but  is  often  not  necessary  since  every  solution  to 
the  FOPDE  is  obtained  by  substituting  all  possible  values  for  a.)  In  the  next  section 
we  will  show  that  with  the  help  of  the  characteristic  equations,  a  complete  integral  can 
found. 

The  general  solution  does  not  comprise  all  solutions  to  a  FOPDE.  The  envelope 
of  a  complete  integral,  the  so  called  singular  integral,  is  a  solution  which  cannot  be 
obtained  from  the  general  solution.  The  equation  of  the  singular  integral,  which  does 
not  contain  any  arbitrary  elements,  is  found  by  eliminating  the  parameters  a  and  b 
from  the  equations: 

<P(i,y,a,6)  ==« 

<Pa[x,  y,  a,  6)  =  0  (A.4.4) 

<Pb{x,y,a,b)  =  0. 

In  fact,  we  do  not  have  to  know  a  complete  integral  of  a  FOPDE  in  order  to  find 
the  singular  solution.  Note  that  for  a  complete  integral  0,  F(x,y,0,0x,0y )  vanishes 
identically  for  all  choices  of  the  parameters  a  and  b.  Differentiating  the  FOPDE  with 
respect  to  a  and  b  we  obtain: 

F*0a  +  Fv0xa  +  F,#vo  =  0  (A.4.5) 

F#0b  +  Fp&xb  +  Fq&yb  =  0. 

As  0  is  a  complete  integral,  A  =  0xa0yb  —  0xb0ya  is  not  equal  to  zero.  Furthermore 
0a  and  0b  are  zero  (equations  (A.4.4))  and  therefore  we  may  derive  the  equation  of  the 
singular  integral  by  eliminating  p  and  q  from: 

Fp  =  0  F,  =  0  F  =  0.  (A.4.6) 

Note  that  equations  (A.4.6)  are  derived  from  (A.4.5).  (Remark:  In  this  case  we  do  not 
assume  that  F£  +  F^  ^  0  as  we  did  to  obtain  the  characteristic  equations.) 

If  a  FOPDE  does  not  contain  the  function  z[x,y )  explicitly,  then  no  singular 
solution  can  exist  as  in  this  case  the  complete  integral  is  of  the  form  [COIH62b,  p.24]: 

z  =0{x,y,a)  +  b  (A.4.7) 

and  the  condition  #6  =  0  cannot  be  fulfilled. 


* 
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A.5.  A  Method  for  Finding  the  Complete  Integral 

In  the  previous  sections  we  showed  two  methods  for  finding  the  solutions  to  a  given 
FOPDE.  Here,  we  explain  how  to  find  the  complete  integral  using  the  characteristic 
equations.  Furthermore  we  describe  how  to  determine  a  one-parameter  subset  from  the 
four-parameter  family  of  solutions  to  the  characteristic  equations. 

First  let  us  discuss  a  special  form  of  a  FOPDE  called  PfafPs  equation: 

f(x,y,z)dx  4-  g(x ,y,z)dy  -f-  h(x,y,z)dz  =  0.  (A.5.I) 

In  the  case  where  h  =  0  and  /  and  g  depend  only  on  x  and  y,  this  equation  degenerates 
to  an  ordinary  differential  equation  called  an  exact  differential  equation: 

f(x,  y)dx  -f  g(x,  y)dy  =  0.  (A.5.2) 

The  equation  is  called  total  if  /  and  g  satisfy  the  integrability  condition: 

fy[x,y)  =  gx{x,y).  (A.5.3) 

In  the  case  of  a  total  differential  equation  it  is  easy  to  find  a  solution  to  (A.5.2).  On 
each  simply  connected  region  we  can  determine  a  function  H[x,  y)  such  that: 

dH  .  ,  dH  .  .  , .  . 

—  =  /(x,  y)  and  ~  =  ff(x,  y).  (A.5.4) 

Then: 

dH  =  f(x,  y)dx  -f  g(x,  y)dy  (A.5.5) 

and  the  equation  dH  =  0  are  both  equivalent  to  (A.5.2).  Thus  H(x,y)  =  c,  where  c 
is  a  constant,  is  a  solution  to  (A.5.2)  and  the  function  H  can  be  found  by  iutegrating 
(A.5.4). 

In  the  case  where  a  FOPDE  is  exact  but  not  total,  an  integration  factor  n[x,y) 
can  be  always  introduced  such  that  the  equation: 

Hfdx-\-  ngdy  =  0  (A.5.6) 

is  total,  i.e.,  such  that  (/i/)y  =  {ng)x-  Equivalently,  n[x,y)  has  to  be  a  solution  to  the 
following  FOPDE  which  we  can  solve  using  the  method  of  characteristic  curves: 

A»(/v  —  9x)  4-  —  Hx9  =  0.  (A.5.7) 

Equation  (A. 5. 1)  is  also  easy  to  solve  if  its  left  hand  side  is  a  total  differential  of  a 
function  H{x,y,z)  i.e.,  if: 

dH  .  dH 
9  *  dy  _  dz' 


Ban 1 


(A. 5.8) 
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For  the  previous  equations  to  hold,  i.e.,  for  the  functions  /,  g  and  h  to  be  the  first 
order  partial  derivatives  of  a  function  H,  it  is  necessary  that  the  following  integrability 
conditions  be  satisfied: 


£/=5£  dg  =  dh  dh=df 

dy  dx  dz  dy  dx  dz'  '  ’  ' 

In  a  simply  connected  region  these  integrability  conditions  assure  the  existence  of  a 
function  H(x,y,z)  which  satisfies  (A.5.8)  and  can  be  calculated  as: 


H{x,y,z)  =  \ 

J  fj 


(zo>Voi*o) 


(fdx  gdy  +  hdz)  -f  C 


(A.5.10) 


where  (xo,yo,Zo)  is  a  fixed  point  and  C  is  a  constant.  Clearly  H[x,y,z)  —  c,  where  c 
is  a  constant,  is  a  solution  to  (A.  5.1). 

In  (A.5.9)  is  not  satisfied,  it  is  again  desirable  to  find  an  integration  factor  fi(x,y,z ) 
such  that  the  expression  nfdx  +  ngdy  -(-  fihdz  is  a  total  differential  of  a  function.  In 
contrast  to  the  case  of  Pfaff’s  equation  in  two  variables,  it  is  not  always  possible  to  find 
such  a  factor.  For  such  a  fi  to  exist,  it  is  necessary  that  the  following  equation  holds: 

f{g,  -  hy)  +  g(hx  -  /,)  +  h{fy  -  gx)  =  0.  (A.5.11) 


It  can  also  be  shown  (simple,  but  tedious)  that  in  a  simply  connected  region  the 
previous  equation  is  sufficient  for  (A.5.1)  to  possess  a  one-parameter  family  of  solutions 
H[x,  y,  z )  =  c.  We  demonstrate  now  how  to  construct  such  a  function  H(x,  y,  z). 

First  consider  the  abbreviated  equation: 


f(x,  y,  z)dx  +  g(x,  y,  z)dy  —  0.  (A.5.12) 

This  is  a  PfafPs  equation  in  the  two  variables  x  and  y,  with  rasa  parameter.  Thus  a 
solution  to  it  can  be  always  found  (if  necessary  with  the  help  of  an  integrating  factor 
M  x,y,z)): 

$(x,  y)  =  u(x,  y,z)  =  c  (A.5.13) 

where  c  is  a  constant.  Note  that: 

x/  =  ^  “d  = 

Now  we  define  a  function  S  which  depends  upon  the  three  variables  x,  y  and  z: 


S{x,y,z)  =  \h— 

oz 


(A.5.15) 
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Using  (A.5  13)  we  redefine  S  as  a  function  of  x,u  and  z,  say: 

T(x,u,z)  =  S(x,y,z).  (A.5.16) 

Suppose  that  T  is  independent  of  x.  Then  we  can  find  H(x,y,z)  by  solving  another 
PfafTs  equation  in  the  variables  u  and  z.  In  fact,  it  is  easy  to  see  that  T  is  independent 
of  x.  One  just  has  to  prove  that  %£  =  0  which  we  do  using  equations  (A.5. 14)  and 
(A.5.15): 

A(U-S)=«„  =  |(X/) 

A(X/1  -  S)  =  uv,  =  £(Xg)  (A.5.17) 

A(x,)  =  ...  =  A(x/). 

Equations  (A.5.17)  written  out  in  full  are: 

Sx  =  h\x  -  /X,  +  \{hx  -  fM)  (A.5. 18) 

— Sy  =  g\M  —  h\y  +  X(j,  —  hv)  (A.5. 19) 

0  =  f\y-  g\x  +  Hfy  -  gx )•  (A.5. 20) 

Multiplying  (A.5. 18)  by  g,  (A.5. 19)  by  /,  and  (A.5.20)  by  h  and  then  adding  up  the 

three  equations  using  condition  (A.5.11)  yields: 

gSx  -  fSy  =  0.  (A.5. 21) 

Differentiating  (A.5.16)  with  respect  to  %  and  y  gives: 

Sx  =  Tx-\-  Tuux  (A.5.22) 

Sy  =  Tutly 

By  combining  (A.5. 14),  (A.5. 21)  and  (A.5.22),  we  obtain: 

0  =  gSx  —  fSy  =  gTx  +  9Tuux  —  fTuuy  =  (A.5.23) 

=  gTx  -j-  \fgTu  —  X/jT^  =  gTx. 

Because  g^Owe  may  conclude  that  Tx  =  0.  Thus  we  rewrite  (A.5.16)  as: 

S(x,  y,  z)  =  T(u,  z).  (A.5.24) 

Equation  (A.5.1),  after  multiplication  by  X  and  using  the  expressions  (A.5. 14)  and 
(A.5.15)  becomes: 

\{fdx  -f-  gdy  -f-  hdz )  =  uxdx  uvdy  -f-  (uM  +  T)dz  =  0 


(A.5.25) 


or,  equivalently: 


du  T(u,  z)dz  =■  0.  (A.5.26) 

This  is  again  a  PfafTs  equation  in  two  variables,  which  can  always  be  solved  using  an 
integrating  factor.  Its  solution  is  ip(u,z)  =  c  where  c  is  a  constant.  Therefore  the 
solution  to  (A.5.1)  is: 

H{x,  y,  z)  =  ip (u(x,  y,  z),  z)  =  c  (A.5.27) 

which  is  a  one-parameter  manifold. 

Finally,  we  describe  a  method  for  finding  a  complete  integral  of  a  general  FOPDE 
F(x,y,  z,p,q)  =  0.  The  basic  idea  is  that  the  total  differential  of  a  solution  z(x,y)  to 
the  FOPDE: 

dz  =  pdx  -f  qdy  (A.5.28) 

can  be  interpreted  as  PfafTs  equation  in  the  variables  x,y  and  z.  In  order  to  do  so,  p 
and  q  must  be  expressed  as  functions  of  x,y  and  z.  So  assume  that  two  functions  / 
and  g  can  be  found  such  that: 

P  =  /(*,  V,  *,  a)  and  q  =  g{x,y,z,a)  (A.5.29) 

where  a  is  an  arbitrary  constant.  Then  the  FOPDE  and  (A.5.11)  (with  h  ~  — 1): 

fg,-gf,-fv  +  9*  =  0  (A.5.30) 

are  satisfied.  In  this  case  the  solution  to  (A.5.12)  is  a  one-parameter  manifold,  but  since 
a  parameter  a  is  already  build  into  the  equation,  this  solution  contains  two  parameters 
and  is  the  complete  integral.  Hence  the  problem  of  finding  the  appropriate  functions 
/  and  g  has  to  be  solved. 

Suppose  that  a  function  G{x,y,z,p,q )  exists  such  that  the  following  equation  can 
be  solved  for  p  and  q  (or  equivalently  for  /  and  g): 

F{x,y,z,p,q)  =  0  (A.5.31) 

G{x,  y,  z,p,  q)  —  a. 

For  this  to  be  possible,  the  inequality  FpGq  —  FqGp  ^  0  has  to  hold.  We  have  to 
prove  that  /  and  g  obtained  in  such  a  fashion  satisfy  (A.5.30)  identically  in  the  three 
variables  x,y  and  z.  Differentiating  the  equations  (A.5.31)  with  respect  to  x,y  and  z 
gives: 

Fx  +  pxFj,  +  g%Fq  =  0  Gx  +  PxGp  +  qxGq  =  0 

Fy  +  PyFp  +  qyFq  =  0  Gy pvGp qyGq  =  0  (A.5.32) 

Ft  +  PtFp  +  qtFq  —  0  Gt PxGp qtGq  =  0. 

After  expressing  p,,  qlt  pv  and  qz  from  the  previous  equations  in  terms  of  the  derivatives 

of  F  and  G  and  substituting  these  expressions  into  (A.5.30)  wc  obtain  a  linear  FOPDE 
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for  the  function  G: 


FPGX  +  FqGy  +  ( pFp  +  <}Fq)Ga  —  ( Fx  +  pFM)Gp  —  (Fy  +  qFt)Gq  =  0.  (A.5.33) 

We  can  solve  this  equation  by  the  method  of  characteristic  curves.  The  appropriate 
system  of  characteristic  equations  is  the  same  as  the  one  for  F{x,y,  z,p,q)  —  0: 

dx:dy:dz:dp:dq  =  Fp:Fq:(pFp  -f-  qFq):  —  {pFM  -f-  Fx):  —  (qFa  Fy).  (A.5.34) 

Only  one  solution  to  these  equations  which  is  independent  of  F  and  contains  at  least 
one  of  the  variables  p  and  q  is  needed.  Such  an  integral  is  the  desired  function  G 
and  will  always  exist  since  the  solution  to  the  characteristic  equations  comprises  a 
four-parameter  family: 

Vi(x,y,z,p,q)  =  Ci  t  =  l,2,3,4.  (A.5.35) 

The  Vi  are  independent  and  at  least  one  of  them  must  contain  either  p  or  q. 

The  method  we  have  described  is  due  to  Lagrange  and  Charpit.  It  has  the 
advantage  over  the  method  of  characteristic  curves  discussed  in  section  A.4  in  that 
one  needs  only  to  find  a  single  integral  of  (A.5.34)  rather  then  a  four-parameter  family 
of  curves. 

A. 6.  The  Initial  Value  Problem  for  Quasi-Linear  FOPDE’s 

In  this  section  we  attack  the  problem  of  how  to  determine  a  particular  solution  of 
a  FOPDE,  once  the  general  solution  is  known.  First  we  consider  a  quasi-linear  FOPDE: 

a(x>  y,  z)p  +  Kx>  y>  z)<i  —  cix<  y>  *)•  (A.6.1) 

In  particular  we  discuss  the  problem  of  how  to  find  an  integral  surface  z(x,  y)  of  (A.6.1) 
which  contains  a  given  curve  C  in  space.  In  the  literature,  e.g.,  [COHI62b,  p.  40],  this 
problem  is  referred  to  as  Cauchy’s  problem.  Clearly  the  following  questions  have  to  be 
answered: 

1)  What  conditions  on  C  are  necessary  to  make  this  problem  is  solvable? 

2)  When  is  such  a  solution  unique? 

Let  C  be  given  by  three  continuous  differentiable  functions  of  a  parameter  t:  x(t),  y{t),  z(t) 
Furthermore,  we  assume  that  the  projection  of  C  onto  the  x-y  plane  (referred  to  as  Co) 
does  not  contain  double  points  and  that  zt2  -f-  y\  0.  If  Co  contains  double  points,  an 
integral  surface  with  self-intersections  is  obtained,  hence,  z  is  not  everywhere  a  single- 
valued  function  of  x  and  y  which  implies  that  along  the  line  of  intersection,  p  and  q 
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are  discontinuous.  Now  to  find  a  solution  to  the  FOPDE  containing  C,  we  construct 
a  characteristic  curve  through  each  point  of  C.  The  equations  for  the  characteristic 
curves  depend  upon  two  parameters: 

x  =  x(a,  t)  y  =  y(a,t)  z  =  z(a,  t).  (A.6.2) 


Note  that  the  functions  x,y,  z  are  still  continuously  differentiable.  To  get  the  equation 
of  the  integral  surface,  we  must  eliminate  the  parameters  s  and  t  from  the  previous 
equations,  i.e.,  we  must  express  s  and  t  in  terms  of  x  and  y.  A  sufficient  condition 
for  this  to  be  possible  is  that  the  functional  determinant  A,  which  is  specified  by  the 
following  equation,  does  not  vanish  along  the  curve  C: 


dx  dy  dy  dx 
da  dt  da  dt 


(A.6.3) 


Using  the  characteristic  equations,  we  rewrite  (A.6.3)  as  A  ==  .  Thus  if 

A  0,  we  may  express  2  as  a  function  of  x  and  y  and  be  assured  that  C  lies  on  the 
surface.  This  solution  is  also  unique  which  is  a  consequence  of  the  following  lemma: 


Lemma  [COHI62b,  p.64]:  Each  characteristic  curve  which  has  one  point  in 
common  with  an  integral  surface  lies  completely  on  this  surface. 


This  lemma  follows  from  the  uniqueness  theorem  for  solutions  to  ordinary  differential 
equations. 

We  can  interpret  the  determinant  A  as  the  outer  product  of  the  two  vectors: 

MI)  “d  MI)  (a-64) 

which  are  the  projections  of  the  tangent  and  the  characteristic  direction  onto  the  x-y 
plane,  respectively.  In  the  special  case  where  A  vanishes  along  C,  these  two  directions 
coincide  and  we  may  deduce  that  C  has  one  of  the  following  three  properties: 

1)  C  is  a  characteristic  curve. 

2)  C  is  the  envelope  of  the  characteristic  curves  (called  an  edge  of  regression). 

3)  Co  is  the  envelope  of  the  projections  of  the  characteristic  curves  onto  the 
x-y  plane. 

We  discuss  when  case  1  occurs  first.  From  A  =  0  we  infer  that: 

1  dx  1  dy 
a  dt  b  dt 


(A.6.5) 
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Using  x (t)  and  y(t)  (from  the  equation  for  C)  in  z( x,  y)  the  following  equation  has  to 
hold  along  C: 

dz  dx  ,  dy  ,  ,  ,  A 

3 f  =  ‘‘*+2»3?  =  a!,  +  6?  =  c'  (A-6'6) 

This  means  that  C  satisfies  the  characteristic  equations  and  is  therefore  a  characteristic 
curve.  In  this  case,  there  exist  infinitely  many  surfaces  through  C  which  satisfy  the 
FOPDE.  To  see  this,  just  choose  another  curve  C  which  has  a  point  P  in  common  with 
C.  Now  to  construct  the  solutions  through  C,  a  characteristic  curve  is  passed  through 
every  point  of  C,  in  particular  through  P.  The  characteristic  curve  through  P  is  C,  thus 
an  integral  surface  through  C  contains  C.  In  this  manner  we  can  construct  infinitely 
many  integral  surfaces  which  contain  C.  They  all  meet  along  the  characteristic  curve 
C  which  therefore  can  be  viewed  as  a  branch  curve. 

One  assumption  made  throughout  should  be  stressed  again  here:  we  require  the 
solutions  to  a  FOPDE  to  be  continuous  and  continuously  differentiable  in  some  neigh¬ 
borhood  of  C.  It  might  be  possible  to  find  a  solution  z  through  C,  along  which  A 
vanishes,  without  C  being  a  characteristic  curve  as  occurs  in  cases  2  or  3  mentioned 
above.  However,  the  derivatives  of  z  are  then  not  continuous  along  C.  This  fact  is 
illustrated  in  the  following  example. 


Example: 

We  wish  to  solve  the  following  equation: 

F{x,  y,  z,  p,  q )  ==  3 {z  —  y)2p  —  q  —  0.  (A.6.7) 

The  characteristic  equations  are: 

3(*  -  y)2 

—1  (A.6.8) 

0. 

The  solution  to  these  equations,  with  the  initial  values  xoi  yo>2o>  ig: 

x  =  (»o  —  yo  +  s)3  +  *o  —  (*o  —  yo)3 

y  =  — « +  y  o  (a.6.9) 

z  —  z0. 

Now  we  impose  the  constraint  that  the  integral  surface  passes  through  C  which  is  given 
by: 

i  =  0  y  =  t  z  =  t.  (A. 6.10) 


dx 

ds 

dy 

ds 

dz 

ds 
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Note  that  C  is  not  a  characteristic  curve.  Setting  the  initial  values  xo,yoi  zo  (i-e->  xtV>z 
for  s  =  0)  to: 

*0  =  0  y0  =  t  zq  —  t,  (A.6.11) 

(A.6.9)  becomes: 

x  =  a3 

y  =  —s  +  t  (A.6.12) 

z  =  t. 


In  this  case  the  determinant  A  is: 

^ ^ _ 3s2 
ds  dt  dt  ds 

Thus  along  the  curve  C  (i.e.,  s  =  0)  A  ==  0. 

There  is,  however,  a  solution  to  the  FOPDE  which  contains  C : 


z  =  x*  +  y. 

Note  that  p  =  Jx— ^  does  not  exist  along  C  (as  x  =  0  there). 


(A.6.13) 


(A.6.14) 


In  the  case  of  a  linear  FOPDE  we  can  make  some  further  statements  about  the 
solution  in  the  case  where  A  vanishes  along  C.  Here,  the  integral  surfaces  are  cylindrical 
surfaces  perpendicular  to  the  x-y  plane,  i.e.,  the  function  defining  an  integral  surface 
is  independent  of  z.  The  linear  FOPDE  is: 

<*(*>  y)p  +  K1*  y)<?  =  c(x>  y)-  (A.6.15) 

Recall  that  A  is  defined  as: 

A  =  Xjj/t  ~  xty9.  (A.6.16) 

Using  the  characteristic  equations: 

x„  =  a 
y»  —  b 

the  following  equation  is  obtained  for  Aa: 

Aa  =  aayt  +  ayat  —  baxt  —  bxat  = 

=  aayt abt  —  baxt  —  atb.  (A.6.18) 

Note  that  if  As  (the  first  order  partial  derivative  of  A  with  respect  to  s)  and  A  vanish 
along  C,  then  A  vanishes  everywhere.  The  proof  of  this  last  assertion  follows  from  the 


(A.6.17) 
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existence  and  uniqueness  theorem  for  ordinary  differential  equations.  By  differentiating 
a  and  b  with  respect  to  s  and  t  and  using  relations  (A.6.17),  equations  (A.6.19)  are 
obtained: 


=  &xa  4" 

at  =  ax xt  -f-  av! It  (A.6.19) 

bt  =  bxa  +  byb 
bt  =  bxxt  +  byyt 

and  so: 

A,  =  (a,  +  bv)  A.  (A.6.20) 

To  express  x  as  a  function  of  y  and  z,  say,  x  —  f(y,z),  we  have  to  assume  that 
y3zt  ~  ZsVt  7^  0  along  C. 

Now,  in  order  to  prove  that  the  integral  surface  2  is  a  cylindrical  surface,  it  suffices 
to  show  that  fz  —  0.  Differentiating  x  with  respect  to  «  and  t  we  derive  the  following 
equations: 


xa  —  fyVa  4”  f aZa 


=  fyVt  4-  f**t- 

(A.6.21) 

Then: 

A  =  f,{z3yt  —  zty3) 

(A.6.22) 

from  which  it  follows  that: 

o 

11 

< 

(A.6.23) 

A. 7.  The  Initial  Value  Problem  for  General  FOPDE’s 


Here  we  pose  the  question:  What  are  the  constraints  necessary  to  determine  a 
solution  to  a  general  FOPDE  uniquely?  Clearly  more  information  than  in  the  quasi- 
linear  case  is  needed  as  now  the  solutions  to  the  characteristic  equations  form  a  three- 
parameter  family  of  curves.  So  let  C  be  a  curve  given  by  i(<),  y[t)  and  z{t )  such  that 
neither  C  nor  its  projection  onto  the  x-y  plane  have  double  points.  Furthermore  p{t) 
and  q(t)  along  C  have  to  be  specified  such  that  the  condition: 


dz  dx  dy 

dt  ~  P~dt  Q~dt 


(A.7.1) 


holds  and  the  FOPDE  (A. 3.1)  (i.e.,  F  =  0)  is  identically  satisfied  in  t.  Thus  the 
functions  x{t),y(t),  z(t),p[t)  and  q{t)  define  an  initial  strip  which  wc  denote  by  Ci. 
From  now  on,  the  procedure  is  very  similar  to  the  one  used  in  solving  the  initial  value 
problem  for  a  quasi-linear  FOPDE.  Through  every  element  of  C\  a  characteristic  strip 
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is  constructed  which  can  be  written  as  x(s,  t),  y(s,  t),  z(s,  t),  p(s,  t)  and  g(s,  <).  To  be  able 
to  express  the  parameters  s  and  t  in  terms  of  x  and  y: 

A  =  x3yt  —  3*ya  =  Fpyt  —  Fqxt  (A.7.2) 

cannot  vanish  identically  along  the  initial  strip.  Assuming  this  holds,  z,p  and  q  can  be 
expressed  in  terms  of  x  and  y.  It  then  remains  to  check  whether  p  and  q  written  in 
such  a  fashion  are  the  first  order  partial  derivatives  of  the  integral  surface  z(x,y).  This 
involves  showing  that  the  quantities  U  and  V: 

U  ~  zt—pxt  —  qyt  (A.7.3) 

V  =  za—pxa  —  qya 

vanish  identically.  Since  we  assumed  that  A  7^  0  we  deduce  from  the  previous  equations 
and: 


0  =  zt  —  zxxt  —  zvyt 
0  =  za  —  zxxa  —  zyya 

that  zx  =  p  and  zy  —  q.  Recall  now  the  characteristic  equations: 


dx  _ 

*  = 


Using  the  first  two  of  these  equations  in  the  last  one,  we  obtain: 

dz  dx  dy 

Ts=pT8+qds 


(A.7.4) 


(A.7.5) 


(A.7.6) 


which  implies  that  V  vanishes  identically.  Now  to  prove  that  also  U  vanishes  identically, 
note  that: 


%at  Paxt  Pxat  Q»yt  QMat 


dU 

ds 
dV 

—  xat  Ptxa  Pxat  QtVa  QVat- 


Subtracting  (A.7.7)  from  (A.7.8)  yields: 


dU  dV 

~0s  ~~  ~dt  =  ~PaXt  ~  PiXi  +  q*yt  ~  q*y *■ 


(A.7.7) 

(A.7.8) 

(A.7.9) 


Taking  into  account  the  characteristic  equations  and  the  fact  that  V  =  0  implies 
Off  =  0,  we  rewrite  the  previous  equation  as: 

dU 

-qJ  —  PtFp  +  qtFq  -)-  xtFx  -(-  VtFy  +  (pit  +  Wt)^*-  (A.7.10) 
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Hpwever,  the  FOPDE  F  —  0  holds  identically  in  s  and  t.  Differentiating  F  with  respect 
to  t  gives: 

Fx it  +  Fyyt  +  Ftzt  -f-  FpPt  +  Fq qt  —  0.  (A.7.11) 

Using  the  previous  equation  in  (A.7.10)  we  obtain: 

d-~  =  — FmU  .  (A.7.12) 


For  any  fixed  t,  (A.7.12)  is  an  ordinary  differential  equation  for  U  as  a  function  of  s 
with  the  solution: 

U{s)  =  U{0)eJo  ~F‘U.  (A.7.13) 

Since  by  assumption,  t/(0)  is  zero,  U  vanishes  everywhere. 

To  summarize  the  previous  results:  given  a  curve  x(t),y(t)  and  z(t)  along  which 
p(t)  and  q(t)  are  known  such  that: 


dz  dx  dy 

dt  ~Pdi+Qdi 
F{x[t),  y(t),  z(t),  p{t),  q(t))  —  0 


(A.7.14) 


A  =  Fpyt  —  Fqxt  7^  0, 


there  exists  a  unique  integral  surface  through  the  initial  strip.  We  obtain  a  unique 
surface  because  the  solution  to  the  characteristic  equations  is  uniquely  determined  by 
their  initial  values. 

The  exceptional  case  where  A  =  0  along  C i  is  analogous  to  the  one  discussed  in  the 
previous  section:  there  are  infinitely  many  integral  surfaces  through  C\  if  and  only  if  it 
is  a  characteristic  strip.  Again  we  can  view  any  characteristic  curve  as  a  branch  element 
since,  on  either  side  of  it,  there  can  be  another  member  of  the  family  of  solutions  to  a 
FOPDE,  while  along  such  a  curve  the  first  derivatives  are  continuous.  Note  that  higher 
order  derivatives  along  the  curve  may  be  discontinuous.  If  C\  is  only  a  focal  strip  along 
which  A  =  0,  then  it  might  be  possible  to  find  an  integral  surface  z  containing  it.  As 
in  the  quasi-linear  case,  this  surface  does  not  have  continuous  derivatives. 

Finally,  suppose  C  degenerates  to  a  point  P  with  coordinates  ( xo,yo,zo ).  Then 
the  strip  condition  is  identically  satisfied  for  all  po  and  qo  which  satisfy  the  FOPDE, 
i.e.,  for  all  po  and  qo  which  determine  the  feasible  tangent  planes  at  P.  So  po  and  qo 
can  be  written  as  functions  of  a  parameter  t.  If  the  quantities  Xo,  Vo,  ^o>  Po(0  and  qo{t) 
are  used  as  initial  values  when  solving  the  characteristic  equations,  a  unique  integral 
surface  which  has  a  conical  singularity  at  P  is  obtained.  It  is  called  the  integral  conoid 
of  the  partial  differential  equation  at  P. 
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Appendix  II 


Transformation 


In  this  appendix  we  show  that  any  solution  to  an  image  irradiance  equation  of  the 
form  (B.l)  can  be  obtained  from  the  solutions  to  an  image  irradiance  equation  of  the 
form  (B.5).  Let 

f{Ap2  +  2 Bpq  +  Cq2  -f  2 Dp  +  2 Eq)  =  E(x,  y )  (B.l) 

be  an  image  irradiance  equation  where  /  is  a  bijection  and  A,  B,  C,  D  and  E  are  real 
constants  such  that  6  >  0  and  AS  <  0,  where  6,  A  and  S  are  defined  by: 


'6=  AC  —  B2 


A  = 


A 

B 

D 


S  =A  + 


B 

C 

E 

C. 


D\ 

E 

ol 


(B.2) 


As  /  is  a  bijection,  equation  (B.l)  and  the  following  transformed  equation  have  the 
same  solutions: 

Ap2  +  2 Bpq  +  Cq 2  +  2Dp  -f  2Eq  =  /-1  (E( x,  y)).  (B.3) 

Equivalently,  we  can  write  the  previous  equation  as: 

Ap 2  +  2 Bpq  -j-  Cq 2  -j-  2 Dp  -)-  2 Eq  =  E{ x,  y). 


(B1) 
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Now  let  z  —  z(x,y)  be  a  solution  to: 

P2  +  92  =  E(x,y)  +  c 
where  x,  y  and  c  are  defined  by: 


B  ,  VAC  -  B 2 

x  =  ~x  H - — - y 

v/C  \/C 

y  =  \/Cx 

c  =  Aa2  +  2  fla/3  +  C/32 


and  a  and  are  defined  by: 


a  = 


CD  — BE 

AC  —  B2 
AE  —  BD 


13  AC-B 2‘ 

Then  z(x,y)  —  z(x[x,y),y(x,y))  —  ax  —  /3y  is  a  solution  to  (B.3). 


(B.5) 


(B.6) 


(B.7) 


Proof:  The  proof  proceeds  by  diagonalization  of  a  quadratic  form.  First  we  express  x 
and  y  as  functions  of  x  and  y: 


Cx-By 

y/C{AC—B2)' 


The  first  order  partial  derivatives  of  z  are  abbreviated  by  p  and  q: 


P  = 


dz 

dx 


(B.8) 


(B.9) 


dz 


In  the  following  equations  we  express  p  and  q  in  terms  of  p  and  q  and  vice  versa: 


P 

9 

P 

9 


VC 


Vac  —  b2 
i 

p— •  —  9- 


—  Q 


B 


sfC  \JC(AC  -  B2) 
B{p  +  a)  +  C{q  +  p) 

Jc 

VAC-B\ .  ,  , 

- — - (p  +  q). 

s/C 


-fj 


(B.10) 
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Thus: 

p2  +  q2  =  £[ B\p  +  a)2  +  2BC(p  +  a)(q  4  fi)  + 

C\q  4  /3)2  4  (AC  -  fl2)(p  4  or)2] 

=  Ap2  4  2Aap  +  Aa 2  4  2Bpq  4  2  Baq  4  2B0p  4  2  Ba/3  4 
Cq2  4  2C0q  4  C/32 

=  £(z,  y )  4  c* 

From  this  last  equation  we  deduce  that: 

Ap2  +  2 Bpq  -f  Cq 2  4  2Dp  4  2Eq  =  £(z,  y) 


(B-ll) 


(B.12) 


which  is  the  same  as  (B.4).  | 


